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Abstract. We present a new formalism to study large-scale structure in the universe. The 
result is a hierarchy (which we call the "Helmholtz Hierarchy") of equations describing the 
phase space statistics of cold dark matter (CDM). The hierarchy features a physical ordering 
parameter which interpolates between the Zel'dovich approximation and fully-Hedged grav- 
itational interactions. The results incorporate the effects of stream crossing. We show that 
the Helmholtz hierarchy is self-consistent and obeys causality to all orders. We present an 
interpretation of the hierarchy in terms of effective particle trajectories. 



1 Introduction 



Understanding the evolution of large-scale structure in the universe is an important ingre- 
dient of present-day cosmology. Neglecting the complications arising from the presence of 
baryons, the evolution of large-scale structure is governed entirely by the dynamics of Cold 
Dark Matter (CDM). Even though the Vlasov and Poisson equations specify completely the 
dynamics of CDM at scales much smaller than the horizon, structure formation is still not 
well-understood since the dynamics is both nonlinear and non-local due to the long range of 
the gravitational interaction. 

CDM evolution in the weakly nonlinear and nonlinear regimes so far has been studied 
efficiently only through numerical simulations. Accurate predictions of the CDM correlation 
functions in the mildly nonlinear regime {k ~ O.lMpc"^ at z 0) will be important for 
comparison with ongoing and future experiments, such as BOSS^ and WFIRST^, targeting 
the baryon acoustic oscillations (BAO) in the matter power spectrum in order to put con- 
straints on dark energy models. The numerical simulations conducted so far trying to assess 
the nonlinear effects on the predictions from linear perturbation theory, have suffered from 
one or more of the following: 1) low mass resolution (thus requiring phenomenological halo 
bias to be used for modeling the galaxy distribution, e.g. [1]); 2) sampling variance (e.g. [2]); 
3) studies with small simulation boxes miss effects from the nonlinear evolution of longer 
wavelength modes [3] ; 4) in numerical simulations one usually fixes all cosmological parame- 
ters to infinite precision, and then varies only the dark energy equation of state, which is not 
adequate for getting a complete handle on the uncertainties (e.g. [2]). 

Therefore, apart from the numerical work, numerous analytical methods have been 
employed to study the CDM power spectrum. However, all analytical expansion schemes 
used so far fail to achieve 1% accuracy for the density power spectrum even in the weakly 
nonlinear regime (see below). This is due to the fact that the convergence properties of 
these expansion schemes rely on the existence of a small parameter, and such a parameter 
simply does not exist in the nonlinear regime. Thus, even for the simplified dynamics in the 
Zel'dovich approximation (ZA) [4] all consistent Eulerian analytical expansion schemes fail 
to recover the density power spectrum when the fractional overdensity becomes close to one 
[5]. 

Moreover, most expansion schemes work in the single-stream approximation. Thus, the 
results of these methods are applicable before shell-crossing, i.e. long before virialization can 
occur. This is done by considering only the first two moments of the Vlasov equation which 
reduce to the usual continuity and Euler equations. Since the higher moments of the one- 
particle distribution function are artificially discarded, one closes the system by introducing 
an equation of state, or equivalently a sound speed, which is set to zero for the CDM (along 
with any anisotropic stress). 

All analytical methods based on the single-stream approximation effectively neglect 
high-A; modes, some of which collapse (and therefore develop shell-crossings) as soon as one 
starts evolving the full Vlasov-Poisson (VP) system. It was recently showed [6] that virialized 
structures decouple completely from large scale modes, just as the collision of two galaxies is 
not affected on galactic scales if all stars were substituted by binaries (i.e. high-A; virialized 
objects). Yet, after averaging over the nonlinear structures, the dynamics of the low-A; modes 
is slightly modified due to the presence of non-virialized structures by introducing a non- 

^ http://www.sdss3.org/surveys/boss.php 
■^http://wfirst. gsfc.nasa.gov 
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vanishing effective sound speed and viscosity [6]. In tlie language of statistical mechanics, 
there is an effective field theory which governs the behavior of the low-A; modes, when one 
integrates out the high-A; degrees of freedom. 

Recently various analytical tricks have been put forward to remedy the convergence 
properties in the weakly nonlinear regime of the standard Eulerian perturbation theory [7]. 
These include renormalized perturbation theory [8], the path-integral approach [9], and the 
renormalization group flow [10]. A critique to these methods, which rely on the single- 
stream approximation, was presented by [11] (see [12] as well), showing that at z = for 
k > O.lMpc"^ shell-crossing may affect the power spectrum at a percent level, which is 
comparable to the projected observational errors of future BAO experiments. A similar result 
was obtained by [13], who demonstrate that all of the existing analytical techniques deviate 
at 1% or more from the correct power spectrum derived from simulations for A; > 0.1 Mpc~^ 
at z = 0. Moreover, [13] find that those methods systematically fail to reproduce the density- 
velocity cross-correlation for the same scales. 

Approximating the CDM as a pressureless perfect fluid implies that one discards any 
short-scale velocity dispersion, which is generated through the nonlinear evolution, as shown 
in e.g. [14]. By obtaining the stress tensor from numerical simulations, it was shown in 
[15] that neglecting the velocity dispersion can lead to 1% effect on the power spectrum for 
k > 0.2Mpc~^ at z = 0. A systematic analysis of different expansion schemes using the ZA 
was performed by [5] who shows that for k > O.lMpc"^ there is an agreement (to better 
than 1%) with the exact nonlinear power only for expansions which are not self-consistent, or 
which add an ansatz for the decay at high k of the density response function, which may not 
be adequate for the gravitational dynamics. Therefore, if one wants to achieve 1% accuracy 
at the scales relevant for BAO measurements, one needs to go beyond the pressureless perfect 
fluid approximation, and consider extending the analysis to include the rich structure of the 
CDM one-particle distribution function in phase space. 

Lagrangian Perturbation Theory (LPT) [16], of which the ZA is a special case (it is 
the lowest order in LPT), tries to do just that: include the phase-space information in the 
calculation of the statistical properties of CDM. However, LPT assumes that the velocity 
field in each CDM stream is irrotational in Eulerian space (see e.g. [17]), which is violated 
after shell crossing. Moreover, it still uses the overdensity as an expansion parameter (albeit 
in Lagrangian space), which becomes 0(1) at the scales relevant to the BAO. Furthermore, 
the "standard" LPT still assumes the single-stream approximation at intermediate steps, 
when deriving the second and higher orders in LPT (e.g. [17]). Thus, alternatives must be 
devised. 

Another machinery to study CDM statistics in phase-space is the BBGKY hierarchy 
[18] which couples the n-point phase-space correlation functions in an infinite hierarchy of 
partial differential equations (PDEs). However, the BBGKY hierarchy suffers from a severe 
closure problem as there is no manifest small physical ordering parameter controlling the 
hierarchy in the mildly non- linear and non- linear regimes (see [18] for further discussion). 

Another line of research was followed by Valageas [19] who used the steepest-descent 
method applied to a large- expansion to obtain equations of motion for the phase-space 
statistics of the CDM. He used an expansion in the linear power spectrum P^, which recovers 
the results of standard perturbation theory, and so does not give insight into the effects of 
stream crossing on the power spectrum.^ 

^We postpone a further discussion of his results to Sections 2.2 and 6.5. 
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The aim of this paper is to start from first principles and obtain equations of motion 
governing the evolution of the phase-space statistics of CDM. This will result in a hierarchy, 
which we call the Helmholtz hierarchy (HH), which is self-consistent to all orders, obeys 
causality, and has a physical ordering parameter. 

One can recast the CDM dynamics into a problem in statistical mechanics, since the 
initial conditions for structure formation are such that the density field (along with the 
particle positions and velocities) is a stochastic random field. Thus, each realization of that 
field can be treated as a microstate and the resulting canonical ensemble can be described 
by a well-defined Gibbs partition function [19]. From the partition function one can obtain 
the Gibbs and Helmholtz free energies, which generate the extended BBGKY hierarchy and 
the Helmholtz hierarchy (both defined below), respectively. 

Since structure formation is an inherently out-of-equilibrium process, no simple relation 
(analogous to the fiuctuation-dissipation theorem) holds between the correlation and response 
functions of the CDM. Instead, the n-point correlation functions and the n-point response 
functions of m-th order are coupled through what we call the "extended" BBGKY hierarchy 
(see Section 3). However, as in the case with the standard BBGKY hierarchy, the extended 
BBGKY hierarchy has no manifest small physical ordering parameter in the mildly non-linear 
and non-linear regimes. 

The Helmholtz hierarchy, similar to the BBGKY hierarchy, describes the phase-space 
statistics of CDM.^ However, we show that the Helmholtz hierarchy is regulated by a phys- 
ical ordering parameter, which is schematically given by the fractional difference between 
the acceleration of test particles as given by a Zel'dovich-type approximation (i.e. the ac- 
celeration is assumed parallel to the velocity at an intermediate moment in time), and their 
corresponding true acceleration due to gravity. Therefore, it effectively interpolates between 
Zel'dovich dynamics and fully-fiedged gravitational dynamics. 

We will see that under a sharp truncation of the HH all n-point correlation functions 
of CDM are generated, in stark contrast to the (extended) BBGKY hierarchy. Combining 
this result with the presence of a physical ordering parameter, we will show that the HH 
ameliorates the closure problem of the (extended) BBGKY hierarchy. 

By constructing the above hierarchy we showed that the effects of stream crossing are 
not non-perturbative as suggested by some authors (e.g. [19]), but are the result of carefully 
treating the initial conditions in phase-space. We show that although the initial overdensities 
follow Gaussian statistics, the initial one-particle distribution function in phase-space has 
highly non-Gaussian features, preserving which is crucial for the consistency of the method, 
and for capturing the correct physics as well. 

One can construct closed-form solutions for the Helmholtz hierarchy at each order, and 
show that they are closely related to a Born-type solution to the Vlasov-Poisson equation. 
By doing that to the second lowest order and reabsorbing some higher order contributions, 
we show that the Helmholtz hierarchy (HH) solutions have a natural interpretation. One 
can think of the solutions in terms of an iterative improvement of the ZA using N-body 
simulations in the following way: 

Initial set-up: One starts with a simulation box with an initially uniform distribution 
of CDM particles. One then allows those particles to move according to the ZA, i.e. on 
exactly solvable straight trajectories. The resulting density field from those particles produces 

*The HH governs the evolution of the functional derivatives of the Helmholtz free energy, which is also 
called the 1 Particle Irreducible (IPI) effective action, F, defined later for the CDM. Those quantities are 
related to the phase-space correlation and response functions of CDM through a Legendre transformation. 
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a gravitational potential ^q, which does not influence the particles, but is stored for the next 
iteration step. 

For n = 1,2, ... do: A new simulation box is constructed which is filled with exactly 
the same particles at the same initial positions as in the initial set-up^. One then solves 
the standard equations of motion for each particle, but for the gravitational acceleration one 
uses^ (— V$n_i). The resulting density field from those particles produces a gravitational 
potential which does not influence the particles, but is stored for the next iteration step. 

At each iteration step, the above scheme produces density fields which gradually inter- 
polate between the ZA and the fully non-linear solution. The gravitational potential resulting 
from the ZA is much smoother than the result in the fully non-linear regime, thus result- 
ing in much smaller accelerations, a, especially at small scales. Those small-scale particle 
accelerations gradually increase with each iteration. Since the usual choice for timesteps in 
N-body simulations goes like l/y^jol, one may hope that even if the integrals involved in 
the Helmholtz solutions are not doable numerically, the above N-body scheme will produce 
the desired solutions using very few iterations with very few timesteps each^ to recover the 
mildly non-linear regime. This can result in a speed-up of N-body simulations targeting that 
regime. The results can have a wider applicability than simply calculating the matter power 
spectrum, since other statistics can be extracted from the simulation boxes in the same way 
one treats the results from fully non-linear N-body simulations. 

Our analysis mainly builds upon the work of [19] and [21]. Thus, most of the technical 
complications in obtaining the HH arise because the action governing the classical dynamics 
of CDM in the presence of the stochastic initial conditions is built up of non-trivial integro- 
differential operators. To aid the reader, we therefore summarize the results at the end of 
each of the more technical sections. 

In Section 2 we start with a review of the Zel'dovich approximation, deriving the phase- 
space response and correlation functions of the CDM in the ZA. We then write down the full 
VP equation and in Section 3 we derive the extended BBGKY hierarchy, which includes the 
response functions of the CDM. Before we can derive the HH, we have to rewrite the Vlasov- 
Poisson system, given the stochastic initial conditions set by inflation, in the form of a path 
integral. This is done in Section 4. Next we review the Non-Perturbative Renormalization 
Group (NPRG) flow equations in Section 5, which we then solve for the CDM to obtain the 
Helmholtz hierarchy in Section 6. We reintroduce the non-Gaussianities of the initial one- 
particle distribution function in Section 7. And in Section 8 we show how one can write the 
HH in terms of an iterative scheme of N-body simulations. We then summarize our results 
in Section 9. 

^This is done so that one works in the same reahzation of the initial conditions. 

®One has a choice of whether to evaluate the acceleration of each particle, (— V$n-i), at the position of 
that particle in the box of the (n — 1)— th iteration, or at the position of that particle in the box of the n— th 
iteration. Further analysis is necessary to establish which approach yields faster converging results. 

^One needs zero timesteps in the case of the ZA, which is exactly solvable. One can also imagine numerous 
other improvements to the above iteration scheme, such as using the second order LPT solution for the 
initial set-up, which has been shown to improve the behavior of full-blown N-body simulations [20]. Another 
possible optimization would be achieved by filtering out the small scales in the initial conditions. Their 
dynamics cannot be captured by the ZA, and would require a large number of iterations to recover using the 
above iterative scheme. Those scales are hardly relevant, because small scale non- linear power should have 
almost no effect on the mildly non-linear scales [18]. 
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2 Preliminaries 



2.1 The Zel'dovich approximation 

We begin our discussion with an overview of the ZA [4]. It is the lowest order approximation 
in Lagrangian perturbation theory and by construction in the Helmholtz Hierarchy as well; 
and it provides the basis for the expansion which is at the root of the HH. In this section 
we derive the one-point and two-point correlation functions in the ZA, as well as the linear 
response function. This will give us a flavor of the phase-space statistics of CDM. Although 
the ZA breaks down in the nonlinear regime, it captures some of the effects of stream crossing 
of CDM, and thus it will give us an insight which is qualitatively different from the results 
in the single-stream approximation. 

In the ZA, the Eulerian comoving coordinates, of a particle are given by 

x{q,r^) = q + D{^)s{q) . (2.1) 

We denote conformal time with ri\ q are the Lagrangian coordinates; s is a stochastic dis- 
placement field, encoding the initial conditions for structure formation; and D is the linear 
growth factor*. The above equation and eq. (2.2) below can be derived in linear theory (e.g. 
[5]). In the ZA they are assumed to hold in the nonlinear regime as well. 
The stochastic displacement field in Fourier space is given by^: 

s{k) = i^5L{k) , (2.2) 

where k is the Fourier wavevector corresponding to x; and 5l = 5{r]j)/D{r]j) is the linear 
fractional matter density perturbation evaluated today, and r]j is some early time in the linear 
regime. By 5 we denote the fractional matter overdensity, 5 = {pM — Pm)/pm-, where pM is 
the average comoving matter density. Throughout this paper we assume that 6l is drawn 
from a Gaussian random field. This assumption can be dropped easily, but it will make the 
analysis much more cumbersome. The linear power spectrum, Pi(/c), of 5l is given by: 

{8L{k)6L{k')) = 6D{k + k')PL{k) , (2.3) 

where 6o{x) is the Dirac delta function; angular brackets denote ensemble averaging. 

We are now in position to write down the one-particle phase space distribution function. 
Before we do that, however, we have to decide on a velocity coordinate. The conjugate 
velocity to x is given by (a denotes the cosmological scale factor) 

dx , 

(2.4) 

such that Sxd^M is the measure in the one-particle phase space^''. However, we would like 
to use a velocity coordinate, such that in the ZA we have an equation of motion given by 



^The linear growth factor is given by the growing solution of d{at))/dri = 4ttGpmD with D{r)o) = 1 today 
(at r]o), where G is the Newton constant, pM is the average comoving matter density, and a is the cosmological 
scale factor. 

®We define the Fourier transform as 

For brevity we use the same notation for quantities in Fourier and in configuration/phase space. 

^"We have dropped a constant factor of the CDM (particle mass)^. And in general, by phase space we mean 
the {x, v) space. 
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dv/dr] = (The reasons for this choice will become apparent when we discuss the full VP 
equation). The coordinate velocity of the CDM particles in the ZA is 

d7] 

where a dot denotes a partial derivative with respect to conformal time. Therefore, we work 
with the rescaled velocity 

v = {aD)~^v, (2.5) 

which for a test particle in the ZA reduces to v = s{q), which is time-independent by 
construction. 

Along with the velocity coordinate, we transform the one-particle distribution itself. 
So, instead of the invariant one-particle phase space distribution function, f{x,v,rj), we use 
f{x,v,rj) defined as: 

f{x,v,r])= (at)) ^ f{x, aDv, 7]) . (2.6) 



Thus, the number of particles in a phase-space element is proportional to 

fd^xd^v = fd^xd^v . 
In (double) Fourier space, this transformation reads: 

= f(.,w = -;t,,) . (2.) 

where w = aDw is the wavevector corresponding to v; while w corresponds to v. 

We can integrate over all streams at a given position x to write the exact one-particle 
phase-space distribution function in the ZA as: 

fix, v,ij) = J d\5D{v - s{q))5D{x -q- D{r,)s{q)) . (2.8) 

Next we Fourier transform / given in (2.8) with respect to both x and v. 

d^xd'^i 



(27r)6 

Thus, the phase-space distribution satisfies the following Vlasov equation: 

f-Dk-d^f = 0. (2.10) 

Now, we are ready to proceed in calculating the correlation and response functions in the 
ZA. 
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2.1.1 Correlation Functions 



Before we proceed to calculate the exact one-point and two-point functions, let us see how / 
relates to the fractional density perturbation, 5, entering in the single-stream approximation. 
We have 



l + 5{x,ri)= d\f{x,v,7]) = d^vf{x,v,rj) = {27rf fix, w = 0,7]) 



(2.11) 



where we used the partially Fourier transformed /. Thus, we find that the CDM density 
power spectrum is given by 



{6{k,7])6{k',7])) = {27rf if {k,w = 0,7]) f{k',w' = 0, rj)) . 



(2.12) 



In order to calculate the above bracket using the expression for /, eq. (2.9), we use that 
for any linear combination, A, of Gaussian random variables, the following equality for the 
ensemble average holds: 



where ()c denotes the cumulants. Thus, we obtain 

(27r)6(/(fc, w, r,) f{k', w', t]')) = doik + k') 



(2.13) 



(2^) 



(2.14) 



xexpj-^ j d^i^^^l{q,i,k,w,w' ,r],r]')^ , where 

l{q, i, k,w,w', 7], r,') = (^ • if + (^ . m)2 + 2 cos(^ ■ q) ■ I) {i ■ m) 

I = w + D{r])k , m = w' — D{r]')k . 

The above expression is the exact two-point function for CDM in the ZA, which includes the 
effects of shell-crossing. Setting w = w' = 0, the above expression reduces to the expression 
for the two-point function of the density perturbations (cf. eq. (105) in [5]). 
The above equation can be integrated in ^ to obtain: 



(271)^ if ik,W, 7]) f{k',w', 7]')) 



5^(fc + fc')e-^('^+™^); 



(27r)3 



(2.15) 



exp 



iqk-lm {x{q) + 7(g)) + 3(i ■ q) (m ■ 9)7(9) 



where 



(^v = ^ dxPiix) , x{q) = ^ /o°° dxPL{x)jo{xq) , j{q) = ^ dxPL{x)j2{xq) 



(2.16) 



where j/ denote the spherical Bessel functions. 

As a consistency check, expanding to linear order in Pl we obtain 



i27r)^ if ik,W, 7]) fik',w',7j')) = 5D{k + k')PL{k) 



{I -k) {m- k) 

1? 



(2.17) 
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which for w = w' = and r] = r]' reduces to {6{k,r])5{k' ,r])) = 6^^\k + k')PL{k)D'^ , which is 
the lowest order expression for the power spectrum in Standard Perturbation Theory (SPT). 
Using (2.9) and (2.13) we can obtain the ensemble averaged one-point function: 



f = TTTTTs- exp 



(2vr) 



(2.18) 



where / = (/). The exponential damping is due to the fact that the large-scale velocity 
dispersion is given by a^. The delta function above encodes the homogeneity and isotropy 
of the universe. 

2.1.2 Linear response function 

The linear-^-^ response function, 6fa/SCb, tells us how = f{ka,Wa,ria) is modified under 
a forcing (^^ = (^{ki„Wh,r]h) at time rjfj. The forcing can be incorporated in the ensemble 
averaged Vlasov equation as follows: 

f-Dk-d^f = C. (2.19) 
Taking the functional derivative of the above equation, we obtain 

^Va-^ - D{r]a)ka • = Soika - kb)6D{Wa - Wi,)5D{Va " %) • (2.20) 

The solution to this equation is 

^ = 0{r,a - r}b)5D{ka - h)6D(^Wa -Wb+ {D{r]a) - D{r^h))ka^ , (2.21) 

where 6{x) is the Heaviside step function. Above we imposed the following conditions [19] 

-j-^ ^ doika - h)6D{wa - Wb) as r]a^r]b + 0'^ 
oCb 

— = for r]b> rja 
oCb 

— = dD{kb)dD{wb) for r]a>rib. (2.22) 

Kb 

The first two equations are causality conditions: information about the forcing C, propagates 
only forward in time. The third equation comes from mass conservation as J d^xd^vf =constant. 

As argued by Valageas [5], the high-/c damping of the linear density response function 
obtained in Renormalized Perturbation Theory (RPT) [3] is due to the erasure of small-scale 
correlations between the initial and final density field by the advection of high-A: structures by 
the random large-scale flows. Since the high-A; decay of the linear density response function 
is due to large-scale bulk motions, it is also present in the ZA [5, 22]. However, the exact 
linear phase-space response function in the ZA, (2.21), does not exhibit such an exponential 
damping at high k. Still, one can recover the exponential damping of the density response 
function by performing an expansion similar to the one done by [16], but applied directly to 
the initial conditions of the phase-space distribution function (see also Section 6.5). 

Having encountered the phase-space statistics of CDM in the simple case of the ZA, we 
are now ready to turn on the full effects of gravity. 



By "linear" we mean that the response function corresponds to the first functional derivative of / in C,. 
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2.2 The Vlasov-Poisson Equation 

In this paper we work with the fuh one-particle distribution function f{x,v,r]) for CDM, as 
defined in (2.6). The evolution of / is governed by the coupled Vlasov and Poisson equations, 
resulting in the VP equation. In later sections, starting from the VP equation, we extract 
the equations of motion for the phase-space statistics of CDM. Before we do that, however, 
in this section we write down the VP equation in a convenient form, suitable for the analysis 
that will follow. 

Since we are working with the velocity defined in (2.5), at early times when the ZA is 
adequate {v ~ s(q) for test particles), we have an equation of motion given by dv/d-q = 0. 
This will result in the expansion required for the construction of the HH. The exact equation 
of motion for a test particle is given by 



where (j) is the Newtonian gravitational potential. The first term above corresponds to the 
Zel'dovich approximation for the acceleration due to gravity at early times, which is parallel 
to the velocity. For modes well inside the Hubble horizon, (f) is given by the Poisson equation: 

k^^ = - — ^(aD\ . (2.24) 
^ aDdr]\ J ^ ^ 

Using (2.4), the Vlasov equation for the invariant one-particle distribution function, f, 
is given by 

df 1 

— = f + -V • a^f - ad^<f> ■d^f = 0. (2.25) 

d7] a 

Transforming (2.25) according to (2.5) and (2.6), and combining with (2.24), in Fourier space 
we obtain 

f - Dk-d^f- (2.26) 



_dln{aD) 
dr] 



/111 ■ k' 
d'k'-^fik', w = 0, i])f{k - k', w, i])-w d^f 



. 



This is the Vlasov-Poisson equation. The parameter e = 1 for the full CDM dynamics, and 
therefore we drop it altogether through most of our discussion. If one sets e = 0, one recovers 
the ZA, (2.10). Up to a change of variables, the above equation is equivalent to eq. (7) in 
[19]. 

Note that "Jeans' swindle" requires that in the above expression we impose [19]: 

Thus, in a homogeneous universe, for which f{k,w,7]) = (27r)~^(5£)(fc), the force term van- 
ishes. 

It will prove extremely useful to use deWitt's notation (without taking into account 
index placement), where subscripts run over all possible labels: i.e. wavevectors, time and 
possibly fields (several of which appear later in the paper); and repeated subscripts imply 



- 10 - 



summation over discrete labels (e.g. field labels), and integration over continuous labels. The 
VP equation is quadratic in = f{ka,Wa,r]a), the quadratic piece coming from the term 
9x4' • dvf where cp is the gravitational potential given by the Poisson equation (2.24). 
Thus, the VP equation, (2.26), can be written as: 

Labfb = Kabcfbfc , (2.28) 

where L and K are integro-differential operators, defined as follows (no summation below): 

Lab = Da[5ab] , where (2.29) 

Da . 5,. - Di,a)ka ■ 5.. + dHairia)Di,a)) . 

5ab = Snika - kb)5D{Wa - Wb)6D{l]a " %) • 

We have defined the operator, Kabc, to be symmetric in its last two indices (no summation 
below) : 



1 dln[a{7]a)D{7]a) 

Kabc = -(27r)^ ^- -D~\r]a) x (2.30) 

2 dr]a 



-5D{Wb)5D{Wa - Wc)5D{ka - fcfe - kc)5D{r]a " 11c)6D{Va " %) + {b ^ c) 



kb ■ w 

Note that the operators above satisfy the following relations: 

Lab = L_a-b, and Kabc = K_a-b-c , (2-31) 

where a negative subscript implies that all Fourier wavevectors with that subscript must 
flip sign. For example, for 6ab = Soika - kb)6D{wa - Wb)6D{r]a - lib) we have 6a -b = 
Soika + kb)5D{'Wa + Wb)5D{'na — Vb)- Note that we use a comma in subscripts just as a 
separator, and not to denote derivatives. 

Note that the term in quadratic brackets in eq. (2.26) vanishes in the ZA, so that 
eq. (2.26) reduces to the VP equation in the ZA, eq. (2.10). Therefore, we can define an 
ordering parameter given by (no summation below): 



dva In a{r]a)D(j]a) M^a • d^^^fa - Kabcfbfc 

Ea = . (2.32) 

max{/a, D(r]a)ka ■ d^^fa} 

This is exactly the physical ordering parameter of the HH. Let us use the ZA expression for 
/, (2.9), and see how e scales. The numerator vanishes at first order in 6. To second order 
in 6 in the ZA, we have: 

numerator of e ~ 5^ In ( at) ) Dwks^ . (2.33) 



The denominator to first order in 6 is given by ~ -Dfc • s. 

In order to proceed we need an estimate of w. The two point function in the ZA can 
be written as the initial conditions propagated by two response functions until the time of 
interest (e.g. (6.44)). Since we are interested in quantities such as the density or velocity 
power spectrum, G is evaluated at w around zero. Then from the response function in the 
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ZA, (2.21), we can see that the largest w we are interested in is^^ w ~ Dk. Therefore, we 
obtain 

e ■ s DSl ^ 6 in the ZA. (2.34) 

In the second equahty above we used eq. (2.2). 

Once we go outside the hnear regime, the estimate above breaks down. However, we 
can still assign a physical significance to e, starting from the fact that 

denominator of e ~ Kf^ . (2.35) 

One can check (2.35) both for virialized objects and for linear perturbations. For linear 
perturbations we have 

\Kf\^{h/a)\ws\r^\f\ , 

where we used w ~ Dk, and D/D ~ a/ a. While for virialized objects, we have Kf^ ~ 
—Dk ■ dwf . An easy way to see that is to realize that for a quasistatic halo / vanishes, and 
the term containing w -dwf is suppressed, as it is proportional to acceleration in the ZA (i.e. 
it measures effects of the Hubble flow on the scale of the halo), while K is proportional to 
the acceleration of gravity in the halo.^'^ 

Combining (2.32) with (2.35) we can see that e can be schematically represented by 
the fractional difference between the acceleration of test particles as given by a ZA (as can 
be seen in the first term in the numerator), and their corresponding true acceleration due 
to gravity (present in the second term in the numerator). The ordering parameter for the 
HH effectively interpolates between Zel'dovich-type dynamics and fully-fiedged gravitational 
interactions. Since K is first order in e, we can see that a truncation to C'(e"'') requires a 
truncation to order O ((X/^)"^), with uk = ng. Note that e in (2.26) was inserted, so that 
we can easily keep track of the order of the expansion: A truncation at 0{i"') is equivalent 
to a truncation at O^e"'), and vice versa. 

It will prove useful to incorporate in the VP equation the initial conditions for / at time 
r]j (no summation below): 

fl,a = f{ka,Wa, 1]l)6D{'na - m) ■ (2.36) 

As we will see in Section 3 the resulting modified VP equation will modify the standard 
BBGKY hierarchy [18] to include the response functions of the CDM dynamics. The modified 
VP equation reads: 

Labfb - Kabcfbfc = fl,a ■ (2.37) 

This way the initial condition for this modified equation is f{r] < rjj) = 0. In the case of 
CDM, f{ka, Wa,r]i) is the phase-space distribution obtained in the linear regime after matter- 
radiation decoupling. Therefore it can be obtained with linear theory from the inflationary 

If for some reason one is interested in the statistics at very large w (e.g. to calculate statistics involving 
/ fv"d?v with high n), the analysis breaks down, since then we would have w ~ max (Z)fe, WqJ- jjj^gj.gg^) ~ 
™of interest ' ^^'^ such a case e may not be sufficiently small for an adequate expansion. 
^•^For the isothermal sphere halo profile, for example, one obtains 

(halo length scale) x a/a ^ 
halo velocity dispersion 
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initial conditions, and thus, the forcing term, //, is a stochastic random variable, which acts 
only at 777. 

The first two cumulants of // are given by: 

/"a = (//,a) and Aab = {fl,afl,b)c , (2.38) 

where angular brackets denote ensemble average over the stochastic initial conditions set by 
inflation, and we denoted the cumulants with ()c. 

Both Ha and A^fe can be evaluated in the ZA by reading off / and (/^)c from eq. (2.15) 
and eq. (2.18), using (2.36) and (2.38). One may think that the appropriate choice for the 
initial conditions is to keep f{i]i) correct only to first order in 6{r]j), as in eq. (2.17), which 
is what [19] uses. However, one should remember that the phase-space initial conditions for 
/ and (/^)c given by (2.18) and (2.14) evaluated at r]j should be kept consistent for w up to 
w ~ Dk (see discussion ca. eq. (2.34)), with D evaluated today (or at the rj of interest) and 
not at r]j. Therefore, expanding eq. (2.15) to first order in Pl is correct to 0{5{r])'^) and not 
to 0{6{rii)'^). Thus, eq. (2.15) breaks down even at weakly non-linear scales, and one should 
use (2.14) instead. 

One should note that it is only the overdensity and not the full one-particle distribution 
function which is a Gaussian random field. The full / is in fact non-Gaussian and keeping 
the non-Gaussian parts will turn out to be crucial for the self-consistency of the formalism. 

Indeed, in the ZA, the (higher n)-point functions are generated at higher order in S{r]j). 
For example, the 3-pt function is generated at 0{6{r]j)^). However, following an analogous 
analysis as above, we can see that the effects from the initial (higher n)-point functions will 
actually kick in not at 0{6{r]j)'^) but at 0{6{r])^) for weakly non-linear and nonlinear scales. 

However, in most of what follows we restrict ourselves to Gaussian initial conditions for 
/ in order to keep the formalism more transparent. Thus for now we will assume that the 
statistics of // are entirely specified by fia and Aab- We will restore the non-Gaussianities in 
// in Section 7. 

Note that Valageas [19] expands Aab to first order in Pl, as in eq. (2.17); and drops the 
non-Gaussianities of //. Therefore, his initial conditions for the 2-pt function of / are only 
correct to first order in w. This means that [19] artificially smooths the initial conditions 
in the momentum direction. The expectation value / Spp^fj oc d"" /{dw^)fi{w = 0) is zero 
for n > 1 under this approximation. Since the higher velocity moments vanish, the higher 
order velocity cumulants must be non-zero. Thus, the approximation in [19] affecting the 
initial conditions, does not reduce to a pressureless fluid with vanishing shear stresses at r]j. 
Because of the structure of the equations, and the fact that Valageas works in a Pl expansion 
in order to obtain a final result for the matter power spectrum, that assumption remains valid 
for all times and is not contained only to the initial conditions. 

3 The extended BBGKY hierarchy and the response functions of the CDM 

In this section we review the derivation of the BBGKY hierarchy [18]. Note that the BBGKY 
hierarchy as written in [18] concerns only equal-time n-point functions. We drop that restric- 
tion in this paper. In the process we will also incorporate the initial conditions in the VP 
equation as done in (2.37), and the result will be what we refer to as the extended BBGKY 
hierarchy for both the correlation and response functions. 
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One can try to solve the VP equation, (2.37), using a Born-type approximation, which 
schematically looks as follows: 

/ = L-'fi + L-'Kff = L-^fj + L-^K{L~^fi){L~^fi) + ■■■ . (3.1) 

Note that we have written the series above so that the operator K takes one argument on 
the left, and two on the right. This series has a straightforward physical interpretation. 
Each term can be represented as a diagram, where several modes of the initial conditions are 
propagated forward in time by the linear operator, L~^. Then, the modes interact one by 
one gravitationally through the 3-vertex, K, and the resulting mode is linearly propagated 
in time. It then interacts with other modes until one obtains the mode of interest (/ at a 
certain moment in time). To obtain the statistics of the phase-space density, one needs to 
multiply one or more /'s, as given by eq. (3.1), and then take the ensemble average. One 
thus obtains that the n-point function is given by a series of integro-differential operators 
acting on the statistics of the initial conditions, //, given in (2.38). This approach quickly 
leads to a mindless proliferation of diagrams. One may try to systematically group different 
diagrams under certain meaningful physical quantities. One way to achieve this is by using 
the extended BBGKY hierarchy, which gives us the equations of motion for the correlation 
and response functions. 

The first equation of the extended BBGKY hierarchy is obtained by taking the ensemble 
average of eq. (2.37): 

Lahfh — Kabcfbfc - KabcGbc = /^a , (3.2) 

where fa = (fa) and Gab = ifafb) - fah = {fafb)c denotes the connected 2-point correlation 
function. For a homogeneous and isotropic universe, when the ensemble averaging is done 
without constraints (such as forcing the origin to be on top of a halo as in [23], for example), 
we can see that the term Kabcfbfc vanishes identically, since it is proportional to dx(j) which 
vanishes for homogeneous / as per (2.27). 

Equation (3.2) is an equation for /„. To find Gab we need to multiply eq. (2.37) by / 
and then take the ensemble average. The resulting equation is 

{Lax '^Kaxyfy)Gxb Kaxy{fxfyfb)c — {fl,afb)c j (3-3) 

where we used eq. (3.2). The three point function is given by multiplying (2.37) by and 
then taking the ensemble average. This yields: 

{Lax '^Kaxyfy){fxfbfc)c '^KaxyGxbGyc Kaxy{fxfyfbfc}c — {fl,afbfc)c ■ (3-4) 

And in general, the n-point function is given by the cumulants up to the (n -|- l)-point 
function. 

Setting // to zero in eq.s (3.2,3.3,3.4) we recover the standard BBGKY hierarchy which 
requires initial conditions to be set for the n-point functions (which is especially straight- 
forward for Gaussian initial conditions). The presence of // in those equations means that 
we have incorporated the initial conditions inside the equations at the cost of introducing 
cumulants of the type {fif^)c which are needed in order to close the hierarchy.^^ Those can 

^*This complication will allow us to transform the extended BBGKY hierarchy to the HH. 
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be expressed via the response functions of the system. We will later show that for the CDM 
with Gaussian initial conditions: 



KxiKx2 ---^Cxr, 



(3.5) 

C=o 



where C is a constant-across-realizations forcing, which corresponds to changing the VP 
equation (2.37) by adding the forcing ( to the stochastic //. Thus, 5'^(/")c/'5C™(C = 0) is 
the n-point response function of m-th order. From now on all derivatives with respect to ( 
are meant to be evaluated at (" = 0, even without stating that explicitly. 

In order to obtain (///)c, which enters in eq. (3.3), we use (3.5). Therefore, after that 
we need an equation for 5f /5C,. Taking the variation with respect to Q of eq. (3.2) we obtain: 

{Lax - '2'Kaxyfy)Tr " Kaxy^^J^ = = ^D{k.a - kb)5D{Wa - Wb)6D{,1la " %) , (3.6) 

where in the last equality we used that adding a non-random to // corresponds identically 
to changing the mean of //. To remind the reader, all derivatives in are to be evaluated at 
^ = 0. Both the above equation and eq. (3.4) depend on 6Gab/SC, which can be obtained by 
varying eq. (3.3). This in turn generates a term 5"^ f /SC,"^ . And so on. 

Finding closure relations for the resulting hierarchy of equations for the correlation and 
response functions is a non-trivial task. One can simply set to zero all n-point correlation 
functions for n above some cutoff. However, in the nonlinear regime this is inconsistent as the 
(higher n)-point correlation functions scale as products of (lower n)-point functions, which 
can be much bigger than 1 [18]. Since nonlinear features are generated at small scales as 
soon as one starts evolving the VP system, such a truncation is unsatisfactory even at early 
times if we want to study the effective CDM dynamics. 

Another option is to use some ansatz for the (higher n)-point functions to close the 
extended BBGKY hierarchy. However, such ansatzes are usually extracted only from obser- 
vations or simulations (e.g. [24]). 

In this paper, we investigate an alternative route. Each contraction entering in the ex- 
tended BBGKY hierarchy between the operators L and and the correlation and response 
functions can be represented by a diagram. We will argue that the HH corresponds to the 
extended BBGKY hierarchy with closure relations equivalent to dropping certain diagrams 
from the expressions for the correlation and response functions. The closure relations ob- 
tained in such a way are equivalent to dropping higher order vertices from the IPI effective 
action. This in turn corresponds to performing infinite resummations of certain diagrams 
entering in the expressions for the correlation and response functions in a way which is easily 
generalizable to arbitrary order in the HH. We refer the reader to eq.s (6.48, 6.49) derived 
below which give an example of a truncation to the extended BBGKY hierarchy obtained by 
a sharp truncation to the Helmholtz hierarchy. 



4 Statistical mechanics formulation of the Vlasov-Poisson equation 

In statistical mechanics any analysis of the statistical properties of a system starts by writing 
down a partition function. In this section we rewrite the expectation value of any functional 
F\f] as a path integral which will allow us to write down the partition function for the CDM. 
We do this by closely following [19]. 
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4.1 Path-integral formulation of the Vlasov-Poisson equation 

Following [19], we apply the Martin-Siggia-Rose (MSR) technique [25] [26] to the stochastic 
differential equation (2.37) giving the evolution of the CDM. The method allows us to rewrite 
the ensemble average of any functional F[f], where / is the solution to a stochastic differential 
equation (in our case (2.37)), into a path-integral form. Such a path-integral formulation 
will allow us to apply techniques borrowed from statistical mechanics to study structure 
formation. 

The ensemble average of a functional, F[f], can be written as: 

(F[/]) = / I?/^F[/]e"5{^.«-/^-)[^'lJ^.''-A''>) , (4.1) 



where any normalization factor is absorbed in D. The above equation is nothing but an aver- 
age of F[f] over the random initial conditions set by inflation, which are taken into account 
by the Gaussian kernel in the integral. The measure is defined as Vfj = Yl^ pdff{x,p) (up 
to a normalization factor), i.e. it goes over the field values on the initial Cauchy surface. 
In the above equation, / is determined by // using the equation of motion (2.37). We can 
rewrite (4.1) as follows 

(F[/]) = I P/,P/det[L,6-i^„,,/e]F[/]e-^(^^--'^")[^"'lo(^^'''-^^) 

^SoiLabfb- Kabcfbfc- fl,a) , (4.2) 

where the delta function imposes the equation of motion. The path-integral measure, T>f, 
goes over all field configurations on all of phase-space and times later than rjj. As shown in 
[19], the determinant reduces to an irrelevant constant, as long as we require that the Green's 
function for 9^^, which enters in Lab, is S^iVa — ilb), where 0"(x) is the Heaviside step function 
with the value at x = given by 6°'{0) = a. As we will show later when we solve the NPRG 
flow equations, this choice for a Green's function will enforce causality^^. 

The next step of the MSR method is to rewrite the delta function as follows 

{F[f]) = j VfiVfVxF[f]e~^^^''''~''''^^^'^^'^b^^''^~''^^~^''^^''''^^~'^'^^''^''^ , (4.3) 

where x is an imaginary auxiliary field. We can now perform the Gaussian integral over // 
and obtain 

{F[f]) = j V<pF[f]e-^^^\ where (4.4) 

SW\ = XaLabfb - XaKabcfbfc " XafJ-a " ^Xa^abXb , (4.5) 

where we combined the fields / and x i^ito the field vector (j) = {f,x)- Thus, the action S[(j)] 
entirely encodes the statistics of the CDM evolution. Equation (4.4) is the final result of this 
section. 



^^Note that 9'^{r]a — T}b) is not invertible, since it is triangular and has zeros on the diagonal r)a = '76- In 
non-equilibrium statistical physics this is dealt with by working with discrete time. What one shows then is 
that actually there are purely imaginary numbers on the diagonal, which can be shown to be irrelevant, (e.g. 
[27]) 
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4.2 Response and correlation functions 

In this section we derive expressions for the response and correlation functions of / based on 
the path-integral, eq. (4.4). The final result will show that the role of the auxiliary field x is 
to generate the response functions of CDM. 

Let us go back to the VP equation, (2.37). We can add a non-random (constant across 
realizations) forcing, ({k,w,r]), on the right hand side of that equation. This modification 
can be taken into account in (4.1) and (4.5) by replacing by /i^ -|- (a- The n-th functional 
derivative of {F[f]) with respect to C gives the response of {F[f]) under a change in the 
forcing (. Using (4.4) with S ^ S — XaCa we find 

6-{F[f]) 



KaiKa2 ■ ■ ■ ^Car, 



= {F[f]Xa,Xa,---XaJ ■ (4.6) 

C=o 

For F[f] = fbifb2 ■ ■ ■ /fem; we see that {F[f]x aiXa2 ' ' ' Xfln) reduces to the 77i-point response 
function of n-th order. 

Causality requires that any effect from the non-random forcing ^ must propagate forward 
in time. This implies that for the bracket above to be non-zero, at least one of the /'s in the 
bracket must follow all ^'s, and hence all x's as well: 

^"(Ai/fe2 ■■■ fbm) ^ 

= (/fei/62 • • • fb^XaiXa2 ---XaJ oc 6*° [ max(% , • • • , r/6„J - max(r/ai , • • • ,VaJ] ■ (4-7) 

As we will later show, the choice for the value of 9(0) above is fixed by the requirement that 
the HH must not violate the above causality condition. 
For F[f] = 1, using {F[f]) = (1) = 1 we obtain 



= {Xa^Xa2 ■ ■ ■ XaJ = , (4.8) 

C=0 



which tells us that the auxiliary field itself is unobservable. From here one can show that 

Sir) 



= irxjc , (4.9) 
C=o 



(shown empirically using Mathematica for a wide range of m and n). 
Let us define the Gibbs partition function (e.g. [28]), Z[j] 

ja4>a\ _ / T)rhp~^\-^\+i 



Z[j] = (e^""^") = J p,/,e-^l'?'l+J'"^- , (4.10) 

where a in ja and (j)a runs over field labels as well. Thus we have ja^a = jf,afa + jx,aXa, 
where jf^a find j^^a are the external sources for / and Xi respectively. Taking functional 
derivatives with respect to ja and then setting ja to zero, we see that as usual Z[j] generates 
the correlation and response functions of the system. 

We are interested in the cumulants of / (i.e. its connected correlation functions), and 
therefore we introduce yet another quantity, the Gibbs free energy^^, given by 

W[j] = In Z[j]. (4.11) 



^^Note that in statistical mechanics the logarithm of the Gibbs partition function (in the Gibbs canonical 
ensemble) gives the Gibbs free energy (up to a proportionality constant); while the Helmholtz free energy is 
proportional to the logarithm of the partition function in the canonical ensemble [28]. In quantum field theory 
(QFT) these definitions are usually swapped (e.g. [29]). We adopt the usual statistical mechanics definition, 
unless otherwise noted. 
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The Gibbs free energy generates the cumulants of 

,<P„---4>.^)c , (4.12) 



— Wai'- 

j=0 



where all subscripts after a semicolon denote functional derivatives as follows: For any func- 
tional A of field z, we have ^[z]; = 6'^A[z]/{6za„ ■ ■ ■ dza^). Note that (i^") is not always 
evaluated at j = 0, as it sometimes denotes the ensemble average in the presence of the 
external source j. It should be clear from the context whether ((/>") is to be evaluated at 
J = or not. 

From (4.8) and (4.12), we conclude that W[j] generates all connected correlation and 
response functions (see eq. (4.9)). As an illustration of the notation, let us consider the 
connected 2-point function for /, finding which is one of the goals of this paper. It can be 

abi 



obtained from Wnh, which can be written in block form as follows^^ 



s 


5 


5 


5 










5 


5 


5 


5 











W,ab\^=^= I ''i'^^'i- "'f'^f- I W[j\ 



{fafb)c 

1^ 



(4.13) 



where we used (4.6) and (4.8). Thus, we converted the problem of finding {fafb)c to a problem 
of finding^^ W-ab- 

As an aside, let us derive eq. (3.5) which gives the cumulant (/"^/p )c entering in the 
extended BBGKY hierarchy discussed in Section 3. In analogy with [5], we express // from 
the VP equation and find 

(r/7) = / Px2^/r (-^ + ^ + A-xye-^['^] , (4.14) 

where (A • x)a = '^axXx- Integrating by parts the above equation n times and taking the 
connected part of the resulting expression, we recover^^ eq. (3.5). Note that eq. (3.5) is exact 
even in the nonlinear regime. As the above derivation shows, this is due entirely to the fact 
that we deal with ensemble averages. In analogy with (3.5) we can express the ensemble 
average of the product of any functional F of f with /J" as 



{F[f]fl,bifl,b2 ■ ■ ■ fl,bm)c = ^biXi^b2X2 • • • Afi 



^^Cxi^Cxz ■■■^Cx 

By analogy, we can calculate 5fi/5C, entering in eq. (3.6): 

Sfl,a 



(4.15) 



C=o 



6Cb 



j VxVfxb + ^« + ^-^-) e""^'"^' = ^ab , (4.16) 



where the last equation is obtained after integration by parts. The above result confirms the 
result from the discussion after eq. (3.6). 



^^For a 2-index operator, such as W-^ab, we use the matrix notation to explicitly show the dependence on the 
field components of (or their respective external sources), i.e. / and X- The (1,1) component in a matrix 
corresponds to the (fafb) component of the operator; (2,2) ~ corresponds to the (XaXb) component; (1,2) - 
to the (faXb) component; and (2,1) - to the (Xafb) component. See for example eq. (4.13). 

^*The matrix W-^ab shown above gives exactly the three non-equilibrium Green's functions which appear in 
the Keldysh formalism for studying non-equilibrium statistical systems (see e.g. [27] for a discussion of the 
relation between the Keldysh formalism and the MSR technique). 

^^We showed that eq. (3.5) can be obtained in the above described way from eq. (4.14) using Mathematica 
for a wide range of m and n. 
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4.3 Summary 

Following [19], in this section we wrote the ensemble average of any functional of /, F[f], in 
a path integral form using the MSR technique applied to the VP equation. We started off by 
writing the ensemble average of F[f] as an integral over the initial phase-space configurations 
of /, weighted by the Gaussian random initial conditions set by inflation. 

The integral can be rewritten as an integral over all phase-space configurations of the 1- 
particle density function of the CDM. We showed that the path integral is over the exponential 
of an action, which is the classical action governing the CDM dynamics in the presence of 
the stochastic initial conditions ~ a result which is standard in the MSR method, and which 
was first derived for the CDM by [19]. We were able to write down that path integral at 
the expense of introducing an auxiliary field, x- This field is unobservable; it generates the 
response functions of the CDM; and by construction it enforces the VP equation of motion 
for / for each realization of the initial conditions. 

Using the path integral formalism, we obtained the Gibbs partition function, which 
generates the correlation functions of the fields, / and x- The logarithm of the partition 
function is the Gibbs free energy, which generates the cumulants of the fields. Therefore, the 
equations of motion for the derivatives of the Gibbs free energy must be exactly given by the 
extended BBGKY hierarchy. 

5 NPRG flow equations 

In this section we will Legendre transform the extended BBGKY hierarchy to obtain the HH. 
We will start by writing down an exact equation for W[j] - the generator of all correlation 
and response functions of the CDM phase-space distribution function. We use the Non- 
Perturbative Renormalization Group (NPRG) fiow based on the effective average action 
method [30] and its variants for obtaining W, which provide an exact non-perturbatively 
correct equation for W [30]. Our discussion of the NPRG flow is based on [31]. 

In the next section we will be able to integrate the NPRG flow equations and will end 
up with the Helmholtz hierarchy. We will argue that the Helmholtz hierarchy is equivalent 
to the extended BBGKY hierarchy discussed in Section 3 after a Legendre transformation. 
Thus, we will argue that a truncation of the NPRG (or equivalently, the Helmholtz) hierarchy 
ameliorates the closure problem of the BBGKY hierarchy. 

5.1 CutofT 

The idea behind the Renormalization Group (RG) flow is to obtain an action which describes 
the "interesting" degrees of freedom (e.g. the low- A: modes), (/><, after one averages over 
("integrates out") the "uninteresting" degrees of freedom (e.g. the high- A; modes), ^>. Let 
us parametrize the boundary between interesting and non-interesting degrees of freedom by 
some cutoff parameter A. The resulting action for the interesting degrees of freedom, Sx, 
can be easily obtained from the partition function Z (4.10) if the functional integration is 
performed only over (py. Thus we obtain 



where we used that (p in (4.10) runs over (/>< and </)>; j< and j> are the external sources 
corresponding to </)< and (/)>, respectively. The flow of Sx with A is contained in Wilson's 




(5.1) 
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approach [32] to renormalization theory. Once ah degrees of freedom are integrated over we 
can see from (4.11) that 5a — > —W. Thus, S\ interpolates between S and —W. This is 
the basic tenet of renormahzation theory: to find a coarse-grained (i.e. high-Zc modes are 
integrated over for a cutoff in k) action Sx governing the behavior of the "interesting" modes, 
starting from the "bare" action, S. 

The transition between integrated ((^> ) and unintegrated ((/>< ) degrees of freedom can be 
made smooth if one introduces "incomplete integration" [32] . To understand how incomplete 
(or "partial" ) integration can be achieved, we give an example with a one-dimensional integral 
over some function f{x): 

fx{x) = I dyf{y)rx{x,y) , with rx{x,y) = e"^("-^)' , (5.2) 

with rx{x,y) called the "cutoff" function (compare the above equation with eq. (5.20)). We 
can see that the "flow" of fx satisfies the following boundary conditions 

/A^oo(a;) = constant x f{x) , and fx^o+ix) = j dyf{y) . (5.3) 

Therefore, fx interpolates between / and its integral. Thus, rx provides exactly the smooth 
transition we need to go from (/)< to </)>. 

Following that prescription, the NPRG flow equations for the CDM dynamics can be 
derived after one modifies the partition function (4.10) by introducing a cutoff function, 7^^^, 
depending on a cutoff parameter. A: 

Zx[3\ = I P0e-^['^]+^''"^"-5<^°^«'>'^'' . (5.4) 



It is important to note that the normalization constant entering in 1)0 above is chosen to 
be A- independent, i.e. it equals 1/Z[0]. We impose this requirement for simplicity, so that 
under differentiation with respect to A, the normalization constant does not produce any 
extra terms. 

Broken into field components, in our case the cutoff can be written as 

n'a, = K I '--^ ° I , (5.5) 

.20 



where TV^ is positive and changes between zero and infinity . Note that (5a, -6 in the (fafb) 



and (XaXb) components of 7^^^ ensures that the cutoff term in Zx[j] has a negative overall 



sign. 

In the exact NPRG flow of Wetterich [30], the cutoff function 7?.^^ is chosen to be a 
momentum^^ cutoff as follows: For X ^ ka the cutoff diverges, TZ'^ — )• oo, while for A <C fca, 
the cutoff is set to zero. However, in what follows we keep the nature of the cutoff 7^^^ and 
the bare action, ^[i;^], completely arbitrary. 



■^"Another requirement for TZa is to rise faster than any power of the fields, so that in the end we are not 
left with incompletely integrated degrees of freedom. However, we do not bother with these details, since in 
the end we take the limit where TZa is a sharp cutoff switching between zero and infinity. 

^^Note that in the context of the RG flow, by "momentum" we mean the wave- vector fc, and not the physical 
momentum of the CDM particles. 
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5.2 Effective action (Helmholtz free energy) 

The NPRG flow equation is easiest to write down for the Helmholtz free energy, T[(p] (defined 
below), which is related to the Gibbs free energy via a Legendre transformation. In the 
absence of a cutoff, the field (p is defined as 

= — F— • (5-6 

Therefore, <p = {(p) is the average ("classical") field evaluated in the presence of an external 
source, j [29]. 

Later (cf. eq. (6.13)) we will see that in the absence of the cutoff, (pa obeys 5T[ip]/6ip = 
j = 0, i.e. one has to extremize F to obtain the dynamics of the classical fields. In quantum 
field theory, finding the extremum of F corresponds to finding the equilibrium field con- 
figuration while including all effects of quantum fiuctuations and at the same time setting 
the external sources to zero. In the context of a classical magnetic system, our restriction 
(j = 0) implies finding the thermodynamic configuration under zero external magnetic field, 
but including the effect of thermal fluctuations. In the context of CDM dynamics, varying 
F corresponds to changes to the system while including all effects of the stochastic initial 
conditions and setting any non-random forcings to zero. The last statement can be seen from 
the fact that any can be absorbed in the initial condition fi, in a similar manner as we did 
for the non-random forcing, (a, discussed ca. eq. (4.6). One can view F as an action, and 
the corresponding equations of motion (cf. eq. (5.10) or eq. (6.13)) are the ones governing 
the dynamics of the classical field, ip, hence the alternative name for F - "effective action". 

In the presence of the cutoff, all thermodynamic quantities acquire a dependence on the 
cutoff parameter, which we put as a subscript explicitly only for Z, W and F for brevity. 
The Gibbs free energy is given by 

Wx = lnZ),, (5.7) 

while Tx is given below. The classical fields and the external sources also depend on A, and 
unless otherwise specified, from now on p and j will denote the A-dependent quantities, fx 
and jA, respectively. We require: 

V^a = {djaWx[ja])\ in the presence of a cutoff 7^^^ , (5.8) 

where a subscript after brackets enclosing a partial derivative indicate a quantity held fixed 
when doing the differentiation. 

It will prove useful (see [31] and our Section 5.4) to modify the usual relation between 
Tx and Wx as follows: 

^Xiy^a] + Wxija] = jafa " \^aKbVb • (5.9) 

As A changes, the value of TZ^ for fixed {ka^Wa^i^a) starts at infinity and then is reduced to 
zero. Thus, F^ and Wx start at some initial conditions in functional space when TZ'^ diverges 
for all {ka,Wa,r]a) which are covered by the path integral in eq. (4.2) (i.e. rja > rjj), and 
"fiow" with A until they reach the true F and W when 7^^ for all (fc^, Wajfja). This flow 
in functional space is the functional RG flow generated by NPRG fiow methods. 
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5.3 Deriving the NPRG flow equation 



Before deriving the NPRG flow equation let us derive several helpful relations between the 
derivatives of the free energies. In the presence of a cutoff, derivatives denoted by a semicolon 
are taken at fixed A, i.e. Tx-a = {dipa^x)x- 

Taking the derivative of (5.9) with respect to (p keeping A fixed we obtain [31] 



where we used 



^X;a = ja - T^abfb , 



d^Wx [j{ip)] ) = {d^j)x{djWx)x = {d^j)x ^ 

X 



(5.10) 



Note that equation (5.10) gives the equation of motion for the classical field (/?. 

The NPRG flow is most easily derived by looking at {dxTx)ip- Taking the derivative of 
(5.9) in A at fixed we obtain 



(dxTx)^ + {dxWx)j + \vJl^ab^b = , 
where we used (5.8) combined with 

{dxWx)^ = {dxWx)j + {dxj)^{djWx)x • 



(5.11) 



(5.12) 



In eq. (5.11), we denoted 7t^^ = dxTZ^^ (Th is is the only exception to the rule that a dot 
represents a partial derivative with respect to conformal time). The derivative {d\W\)j can 
be obtained by working with (5.4) and (5.7). One obtains 



{dxWx)j 



1 



~ 2 

Combining (5.11) and (5.13) we end up with 

{dxTx)^ 



i^abWx-M + 'Pb^a] ■ 



\i^abWx-M 



(5.13) 



(5.14) 



Now let us express Wx-ba using F^. Taking the derivative with respect to ip at constant 
A of (5.10) and using (5.8) we obtain 



5ab = {Tx,ax + 'KxWx^xb = Wx,ax{Vx;xb + 7^i) , and therefore 



Wx, 



(5.15) 



where a superscript {n) denotes the n-th functional derivative; and TZ^ is a shorthand for 
Combining (5.14) with (5.15) we end up with [30]: 



'^xy 



idxTx\ 



.(2) 



(5.16) 



ba 



Equation (5.16) is the final result of this section. It is a NPRG flow equation as derived 
by [30] in the framework of the effective average action method. Taking functional derivatives 
of the NPRG equation (5.16) with respect to the field ip at fixed A one obtains an infinite 
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hierarchy of equations for the functional derivatives F^"^ of r;^. Given F^"^ one can derive 
which is our final goal when evaluated at the final A. 
As an example, the functional flow for F^^^ is given by: 

(dxTx-a)^ = -^Tx-a^y {Wx;^,nLWx;wy} , (5-17) 

where is always to be understood as given by (5.15). Neglecting F^'^^ the flow equation 

(2) 

for Fj^ ' is given by 

X;ab)ip ^ X;awx^ X;byz Wx;u,y [Wx;xJluvWx,vz] + ©(F^^^) . (5.18) 

And in general the flow of F^"^ is determined by all functional derivatives of Tx up to r^"^^^ . 

Note that in principle one can work with the NPRG flow equation written for Wx (5.13) 
which in a full treatment is equivalent to the NPRG flow equation for Vx- However, Wx 
and Tx are related via a nonlinear transformation (5.9). Thus, the functional hierarchies for 
F^"^ and do not result in equivalent dynamics under a truncation at one and the same 
order n. From QFT we know that W generates all connected diagrams, while F generates 
all 1-particle irreducible (IPI) diagrams^^ (i.e. those connected diagrams which cannot be 
disconnected by a removal of a single internal line). Thus, one can expect that a truncation 

(n) (n) 

of the hierarchy will give more accurate results than a truncation of the Wj^ hierarchy 
at the same order in the functional expansion. 

5.4 Initial conditions for the NPRG ^aw 

The initial conditions for the NPRG flow are set at the initial value of A = Aj when 7^^* 
diverges for all {ka, Wa, rja > rji). Using (5.4), (5.7) and (5.9), we can find the initial conditions 
for Fa; by first writing: 

exp (-Fa) = exp Iwx- jafa + ^VaT^ab'Pb I (5.19) 



P(/>exp 



-S[4>] - ]^(t>oKab<i>b + ja4>a - ja^a + ^^faT^abfb 



We can substitute (p ^ <j) -\- in the above equation and combine with (5.10) to obtain 



exp (-Fa) = j 1)0 exp 



-S[ct) + ip]+ 0aFA;a - ^(f'a'R'abfPb 



(5.20) 



When TZ^ diverges for all {ka,Wa,r]a > rjj), exp[— 7^''^(/>^] becomes a delta function, and 
therefore we obtain the initial conditions for the NPRG flow [31]: 

rx,[^a] = S[ipa] (5.21) 

up to an irrelevant additive constant. The above equation holds for r]a = rjj as well. This 
can be seen from the fact that at A = Aj for rja = rfj the cutoff diverges for all rj > ru; and 
the fact that the path integral does not go over r]i and therefore for that rja we have 4> ^ (p 
in (5.19). The above simple result, eq. (5.21), explains the modification that we made in 
eq. (5.9) to the usual definition of F. 



■^^One can extend the above analysis and consider the NPRG flow of TV-particle irreducible efl'ective actions 
(e.g. [33]). However, they generate further complications, and therefore we do not consider them here. 
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5.5 Summary 

In this section we reviewed the formahsm behind the NPRG flow equation (5.16) obtained 
first in [30]. The NPRG flow describes the exact, nonperturbatively correct evolution of the 
effective action, Tx, in functional space. The initial condition for the flow with A is the "bare" 
action, i.e. FaJ^^] = S[ip] (in the case of CDM it is given in eq. (4.5)). When the cutoff A 
reaches a final value. A/, for which the cutoff function, TZ^, vanishes for all {ka,Wa,rja)-, we 
recover the full effective action. Thus, the flow of interpolates between the bare action, 
S, and the full effective action, V. Note that the equations of motion for the functional 
derivatives, T^'^\ are precisely the Helmholtz hierarchy (We postpone further discussion of 
the HH until the next section). 

The NPRG flow can be written as an infinite hierarchy of equations governing the 
evolution of the functional derivatives F]^ of F;^. Truncating the hierarchy at some order n 
in the functional expansion provides us with solutions which are functional approximations 
to the true F^"). Once we have F^") up to some n, we can obtain W^^^ (by differentiating 
eq. (5.15) and using (5.8)) which contains the correlation and response functions of CDM. 

The final result for the full effective action does not depend on TZ^ (as long as it di- 
verges/decays fast enough in its two regimes), but the choice of TZ"^ does affect the quality of 
the approximation once a truncation of the hierarchy is enforced. Thus, certain optimization 
schemes have been developed to deal with the problem [33]. Moreover, certain choices for 
7^'^, such as a sharp cutoff, allow for the analytical integration of the flow. In what follows 
we concentrate on using a sharp cutoff, which will allow us to avoid resorting to numerical 
integrations as much as possible. 

Following the original idea of [32] , the standard choice for a cutoff is a cutoff in momen- 
tum space (see footnote 21). The utility of a momentum cutoff is that for A between Aj and 
A/, it provides us with effective equations of motion for the coarse-grained field (i.e. after the 
high-fe modes have been integrated out) derived from F^. Thus, the high- A; physics influences 
the low-A; physics by modifying the low-Zc equations of motion. In the context of CDM, this 
results in an effective sound speed and viscosity for the long-wavelength cosmological fluid 
[6]; or after introducing certain constraints, this may result in effective Fokker-Planck terms 
in the effective equations of motion (see e.g. [23] for an approximate effective equation for / 
describing an ensemble averaged halo proflle). However, the result of [6] is a result for the 
effective equation of motion for / in a single realization. Extracting from that the nonlinear 
evolution of the statistical properties of CDM involves either a numerical simulation, or an 
analysis similar to the one presented here. Moreover, the parameters of the effective equation 
of motion for CDM still have to be extracted from numerical simulations for example [6]. 
One can try to obtain those parameters from the evolution equations for the correlation and 
response functions obtained by the brute-force application of a NPRG flow. However, under 
a cutoff in k the resulting equations obtained from the NPRG flow are not manifestly causal 
(see footnote 29). That is not a problem when the method is applied to systems in equi- 
librium, since they are invariant under time translations. However, in the non-equilibrium 
setting of CDM dynamics the interpretation of these results becomes severely convoluted. 

Thus, alternative approaches have been developed for non-equilibrium systems based 
on the NPRG flow, with a cutoff function not in momentum but in time [21]. In the context 
of CDM, this method works by allowing the coupling of modes to take place up to a cutoff in 
time. This cutoff is moved from the initial time to the time, which one is interested in. Such 
a choice for the cutoff function generates a manifestly causal NPRG flow. The end result of 
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such a procedure are time evolution equations for F^") — the Helmholtz hierarchy. 

In the next section we concentrate on the NPRG flow with a sharp temporal cutoff 
which will allow us to integrate the NPRG flow, and obtain the HH. 

6 Integrating the NPRG equations for the CDM. The Helmholtz Hierar- 
chy. 

In this section we are going to integrate the NPRG flow equations. We call the resulting 
hierarchy the "Helmholtz hierarchy", since it is a hierarchy for the functional derivatives of 
the Helmholtz free energy, F. As discussed in Section 5.5, unlike standard NPRG methods, we 
are going to use a temporal cutoff as proposed in [21]. The temporal cutoff can be interpreted 
as follows. The nonlinear evolution couples pairs of modes at different times (compare with 
the discussion ca. eq. (3.1)). In the NPRG with the temporal cutoff, the modes are allowed to 
interact only until the cutoff time. As the cutoff is moved forward in time, the gravitational 
interactions are allowed to take place until later times. The discussion in this section closely 
follows [21], since the path integral formulation of the CDM evolution is trivially related to 
the non-equilibrium formalism developed in [21] by a standard change of variables, usually 
referred to as Keldysh rotation. 

We spend a big part of this section to show that the sharp temporal cutoff results in a 
hierarchy which is consistent with causality and the fact that the x field must be unobservable. 
These results must hold in the full nonperturbative treatment, but we perform that analysis 
to show that they also hold for our choice of cutoff after a truncation of the NPRG hierarchy. 

6.1 Causality of the flow I 

Having introduced the NPRG flow equation in Section 5, our starting point is to write down 
our choice for a cutoff function: 

" \0 otherwise 

such that 7^^ oc 0^{r]a — r), which matches the choice made in [21]. Moreover, we adopted 
the notation r for the temporal cutoff from [21] in order to highlight that our analysis from 
now on will be based on the sharp temporal cutoff, and will no longer be valid for arbitrary 
cutoff functions. 

The utility of the choice (6.1) becomes apparent once we consider Wr:aia2---a,„ which is 
a sum of products of ((/>") (ni > n), with the expectation values taken in the presence of a 
fixed j: 

U") = L^^ll^ ^ (6.2) 

- y'ai • • • ^ar J P(^e-^[*1+J-'^'^ - V^ai • • ■ far [(Pa^+i ' ' ' (PaJ ■ 

The last line of eq. (6.2) is written for the case when r of the ip^s are inside the cutoff^^. In 
(6.2) we used (5.10). 



'By ipa being "inside" or "outside" the cutoff we mean that TZ^ — i- oo or 0, respectively. 
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From eq. (6.2) one can see that as long as r > 0, the connected m-point function 
(for m > 1) vanishes, and thus W^"^^ = 0. Therefore, if TZl^^' oo for any of the r/j 
(i = l,...,m), Wr>r]i;aia2-a„^ Vanishes. Since the bare action for the CDM contains only 
equal time operators (see eq. (4.5)), as long as r is larger than max(r/ai,--- ,Vam)^ (0™) 
(m > 0) is T-independent. Thus, we can see that for n = 1: 

'PT>-na,a = ^r='na,a , (6.3) 

and for n > 1: 

Wr>Ti[;aia2 -a„ = Wraj^...a„;a,ia.2---an^ {t — Tai---a„) (6-4) 
= W-aia2-aJ^{r - Taj-aJ , with Ta^...a„ = m.ax{r]a^,- ■ ■ ,rjaj , 

where we have 6^{t — Ta^-a^) = Y\7=i ^^i'^ ~ Va^)- The above equation trivially implies 

{drWr>rjj.aia2-aJj = ;aia2-a„ (t - Ta^-aJ for n > 2 , (6.5) 

where one must be careful when integrating (6.5) to use the correct step function as given in 
(6.4). Thus, the flow of Wr"^ freezes after r exceeds the maximum of the time arguments of 

(n) 

the response and correlation functions entering in W-i as required by causality. 
The flow for F^*^^ must have the same causal property, i.e. 

{drTr-aua2,-,aJ-fi = , for T > Ta^...a„ ■ (6.6) 

Combining (5.15) and (6.4) we find that in order for the above equation to hold we must 
require that 

Wr;axKyWr;yb = -Wr^^-abSnir " Tab) ■ (6.7) 

The above equation will allow us to integrate the NPRG flow. In Section 6.4 we analyze in 
detail the causal structure of the resulting flow equations and their solutions. In the end we 
show that the resulting solutions for Ft"^ are consistent with eq. (6.6), and therefore causality 
is preserved. 

6.2 Solving the NPRG flow equations. The Helmholtz hierarchy. 

Following the results from Section 5.4 the initial conditions for the NPRG flow with a tem- 
poral cutoff are given by 

fH[^.] = 5W[^,] . (6.8) 

Note that above we restored the cutoff index for From (4.5) and (6.8) one can see that 
the initial conditions for the flow are such that the only non-zero functional derivatives 5^"^ 
are for n < 3. 

From (4.5) we find 

^'ni'jT,a ~ XT,xLxa ~ '^XT,xKxyafT,y 

^r]j;XT,a -^axfr^x ^axyfr^xfr^y l^a ^aa;Xr,x ■ (^•9) 
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(2) 

The initial condition, Ff^r^^, obtained from (4.5) is: 

(6.10) 



rjj;ab 



yLab — 2KaxbfT,x ~^ab J 



Later we show (cf. eq. (6.27)) that Xt,x can be consistently set to zero. 

(3) 

Finally, the initial condition, T^ir^j, obtained from (4.5) is T^^.^bc = unless one of 
the 4>i {i = a,b,c) fields is and the other two are /. In that case Trj^-abc = —'^K^^^^^f,-^, 
where 7r(a6c) is the permutation of {ahc) that puts the x field index first. For example, 

^VrJaXbfc = -'^Kn{abc) = -'^Kbac- 

The generation of any non-trivial Fokker-Planck-like terms in the equation for the aver- 
age field, {(/)), is captured by the functional flow at second order in the functional expansion. 
Moreover, any nontrivial evolution of the bare vertex Kabc is captured at third order. Thus, 
the truncation we impose on the NPRG flow is 

r("^4) = rj';^^) = . (6.11) 

The truncation above is simply for convenience. The method can be easily extended to 
include nonzero T^-^\ As P!^"-'*^ involves at least n bare vertices^^, K, this means that 

r(^r.>4) _ ^ ^g_^2) 

including higher order terms, where e is the HH ordering parameter, (2.32). Thus, (6.11) 
corresponds to a truncation at 0{e'^). 

Combining (5.16), (6.4), (6.7), (6.8) and (6.11) it is straightforward to integrate the 
NPRG flow. We only need to specify a boundary condition for P^^^ at Ta- The cutoff 
function vanishes at Ta, and from (5.10), setting the external sources to zero (as per the 
discussion in Section 5.2), we obtain 

rr.;a[¥^r.,a] =0 , (6.13) 

which is the equation of motion for the average field, ip. Thus, from the first order NPRG 
equation, combined with (6.13) we obtain the first equation of the Helmholtz hierarchy: 

^m-Afr^,a] = -\rr,,;al2W;2ie\Ta - T12) , (6.14) 

where we used that Wr^^-ab = W-ab- We denoted the dummy integration subscripts above 
with Arabic numerals, which will make the subscripts below easier to follow. Combining 
equations (6.9) and (6.14) we obtain an equation for (pr^a which is necessary for obtaining 
^rii;ab- We show later (cf. eq. (6.27)) that in order to be consistent with causality we need 
XT,a = 0. 

(2) 

For Pr we obtain 



^Tab;ab — ^rii;ab — — ■^^f :aU^f ;b34W;23W-4lO^ {Tab " t] 



(6.15) 



*This can be seen by considering the number of F*"^^ 's entering in the n-th order of the NPRG flow. 
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We will prove that the term 0{T^'^^) in the above equation vanishes identically due to causality 
and the structure of the bare action. And in general, for n > 1 we will find that F^") depends 
only on the F's up to r'^"'"'"^^ and not on r'^""''^^ We show this in Appendix A. 

(3) 

Finally, for Ft- , the flow equation can be integrated to give the third Helmholtz equation 

^T;abc-^r]i;abc = (6.16) 
1, 



2-rf;al2Ff;fe34rf;c56^;23^;45^;610n^ - t) + (6 O c) 



+ (!?(F(4)) , 

^=''"123456 



with r]! < T < Tabc, the last inequality arising from eq. (6.6). Again, we used that the 
term 0{T^^^) vanishes identically as per Appendix A. The causal structure of the Helmholtz 
hierarchy is discussed in detail in Section 6.4. There we show that the above equations are 
consistent with (6.6), which is the reason why we integrated the flow up to the maximum 
time argument. 

The truncated hierarchy above is not closed unless we specify W-ab- Combining (5.15), 
(6.1) and (6.4) for time arguments larger than r]j we obtain 

W-ab = Wr^^,ab = "^r^.^ax^ r^,;xyW^^^,yb + VF^-^^iax^^y" = ^r^^^ax^ r,^,xyWr^^,yb 

= Wr,^.axTr^y-xyWry,-yb6^{Tab - T^y) ■ (6.17) 

To obtain the third equality above we used that Wr^^^xy oc ijab — Txy) and ni^^ (X e^{T]x - 
Tab)^D{'nx — Vy)-! paying attention to the superscripts of the 6'-functions. For the last equality 
above we used (6.6) and (6.4). Using (6.4) and (6.17), the second equation, eq. (6.15), in the 
hierarchy can be written as 

W,ab = Wr^^,ab = W.axT ^W.ybO^ [Tab " T^y) " (6.18) 
- ]^W.ax'^r-^xl2W.2'iW.iiTr-yMW-yb 9^ {Tab " Txy)9^ {T^y - f) 



Note that the structure of the above equations is the same as that of the flow equations of 
[21], where the sharp temporal cutoff was first proposed. 

Equations (6.14), (6.15), (6.16) and (6.18) are the final result of this subsection. The 
first three of them represent the first three equations of the HH. Showing that the solution of 
(6.18) for Wr^^-ab depends only on the statistical properties of the system at previous times 
will be done after we write down the flow equations in terms of Feynman diagrams. 

6.3 Diagrammatic approach 

In this section we introduce a diagrammatic approach for writing down the NPRG and 
Helmholtz hierarchies. It will allow us to simplify the flow equations; to easily check all 
causality constraints; and to show that the auxiliary field remains unobservable after obtain- 
ing the truncated Helmholtz hierarchy with the choice of a cutoff made in sec. 6.1. 

The NPRG flow equations can be easily cast in Feynman diagrams. Indices correspond- 
ing to X we write as wiggly lines; while those corresponding to / - with straight lines. The 
functional derivatives of W we denote with black lines, while those of F - with gray lines and 
gray dots/blobs. As an example, consider W^"^^ which can be written as 
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a 



\ 



Gr 



ab 



Rriab 



R 



'T;ab 



b b 



b a 



a 



J 



(6.19) 



The above equation combined with (4.13) defines the hnear response function, Rr-ab, and 
the two-point correlation function Gr-ab- The arrow above Rr;ab indicates the direction in 
time of propagation of the linear response (i.e. % > ija in this case); and thus we define 

^T;ba = Rriab- 

To give an example for F^^^ , the vertex T,-.^^ is represented by: 



r - F F = a 

T;Xafbfc 




(6.20) 



Note that and Wr""^ are symmetric, since the functional derivatives commute, e.g. 

r - - - = V - - - 

T;Xafbfc T;fi,Xafc' 

6.4 Causality of the flov^ II 

Now let us show explicitly using the diagrammatic approach introduced in the last section 
that the truncated NPRG/Helmholtz hierarchy is causal and that the auxiliary field is un- 
observable. 

(n) 

To show that, first we prove several properties of Tr . Let us for now assume that all 
IPX vertices whose external legs are all straight lines (corresponding to functional derivatives 
in /) are zero: 



(6.21) 



We will prove this relation shortly. Using (6.17), we see that since T^.j_^j'^ = 0, we obtain 
that the wiggly propagator vanishes: 



(6.22) 



where for brevity we used subscripts x instead of j,^. We will use the same short-hand nota- 
tion for W from now on. Then, we can show that the following relation holds: 




cx (max(r?ci, • • ■ ,Vcn) - max(r?/i, • • •,??/„)) 



(6.23) 



- 29 - 



where z = if the corresponding bare vertex vanishes, i.e. T;^j = 0; and z = 1 otherwise. 
Let us check that the above relation is consistent with the Helmholtz hierarchy. From (6.21) 
we know that there must be at least one wiggly external leg to any IPI diagram. Thus, we 
can see that the IPI m-point function with at least one wiggly and one straight external leg 
is given by 

\ 




where a dashed line stands for either a wiggly or a straight line. The blob represents the 
IPI vertex, r^'"'*. In the above equation we used (6.22) which tells us that to any wiggly 
external leg of a IPI diagram one can attach only the response function S^J^. The dashed 
lines drawn inside the large blob can be attached to either internal propagators, or can be 
external legs. 

The causality condition (4.7) combined with (4.6) implies 

li^^b OC d\7]a - lib) , (6.25) 

which will be checked to be consistent with the NPRG flow further below. According to 

(6.23) , without loss of generality, we can choose each wiggly leg shown inside the blob in 

(6.24) to correspond to the maximum time of the corresponding vertex. For the rightmost 
vertex in the blob we have two choices: the external leg of latest time of that vertex can 
either attach to an internal line in the blob, or can correspond to the external leg with time 
r]f,. In the former case, the chain of ^i^b propagators must eventually close on itself, thus 
forming a closed internal time-loop, which must vanish as per (6.25). In the latter case, using 

(6.25) we can write: 

Vb > Vn+i > Vn > ■ ■ ■ > lu > m > m > m >'na ■ (6.26) 

Thus, the NPRG flow does not violate (6.23) if it holds for some r. Equation (6.23) holds 
for T = rji since the bare action includes only equal-time operators. Therefore, by induction, 
we conclude that (6.23) holds for all r. 

We are now in a position to prove eq. (6.21). The proof goes in complete analogy with 
the above proof of eq. (6.23). We can write the vertex T^.j_^j^,,,j^ as in (6.24), except rjb must 
correspond to a straight external line. Since all external lines must be straight, the chain 
of % ■ab propagators inside the blob must eventually close on itself, which means that the 
diagram vanishes^^, if it vanishes at some r. Since eq. (6.21) holds at r = r/j if Xt = 0, by 
induction it holds for all r. 

^^For this to be the case, note that it is crucial that 9{Q) = in (6.25) 
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So, now we have to see whether Xt vanishes. Its equation of motion is given by the / 
component of eq. (6.14). According to Appendix A, the 3-point IPI vertex entering in (6.14) 
reduces to the bare vertex. Since the bare vertex contains only equal-time operators, from 
(6.21,6.22,6.23,6.25) we see that the right hand side of (6.14) must vanish, and therefore 

Xr = , (6.27) 

as required by the discussion above. 

Next we want to prove that the truncated NPRG hierarchy indeed satisfies the causality 
condition, eq. (4.7), which can be rewritten using (4.6,4.12) as: 

^r;f\,-,hmxa„-,xa^ 0^ [^Mvb„ ' ' ' , ) " max(r/„, , • • • ,r/„J] . (6.28) 

To obtain Wr"''^ from {m < n), one has to take functional derivatives of (5.15) with 

respect to j at fixed r. As an example, for n = 3 we obtain 

Wr;abc = -Wr-axWr-byWr;cz^r;xyz ■ (6.29) 

Since W generates all connected diagrams, while F generates all IPX diagrams, we can see that 
in order to obtain a cumulant or a response function, Wt^\ we must chain one or more 

(2) 

with internal propagators and attach propagators, Wt , to the external legs of the resulting 
diagram. Therefore, if we want to obtain one of the response functions, W^.^, ^, i i , 

' ' ^iXai ■■'Xa; /6i '"/6m ' 

from (6.22) we can see that each of the external rja (corresponding to the derivatives of W 
in Xas) must belong to a propagator h,r-wa with r]a < Vw, where w attaches to the straight 
external leg of a IPX vertex. According to (6.21), that vertex has at least one external wiggly 
leg with time rjx, such that rjx > Vw- That leg must attach to ^t^^ with r]y > 7]x, which is 
either external or in turn attaches to another IPX vertex, thus forming a chain of causally 
ordered IPX vertices. Xn the end, this causally ordered sequence implies that 175^ > r]a^ for 
some r = 1, . . . ,m and all s = 1, . . . , L Thus, we showed that eq. (6.28), and hence eq. (6.25), 
is consistent with the NPRG flow, and by induction it must hold. Thus, we showed that the 
NPRG flow is causal. 

Next we must show that the field Xt is unobservable. We already showed that its mean 
value is zero, eq. (6.27). Now we have to show that all of its cumulants vanish. Combining 
(4.8) and (4.12), that is equivalent to showing that: 

Wr;xu-,Xn=0. (6.30) 

Performing the same analysis as we did above to show the causality of the fiow, we can see 
that all X cumulants must contain a IPX diagram with external legs which are all straight. 
That diagram vanishes, (6.21), and therefore we conclude that indeed eq.(6.30) holds. 

The whole discussion above is self-consistent and consistent with the NPRG flow if both 
the bare action contains equal time operators, and eq. (6.6) holds. The former condition holds 
for the CDM. The latter condition leads to eq. (6.7), which introduces the delta function in 
the NPRG fiow equations, which allowed us to integrate the flow. 

So, we are left to prove eq. (6.6). If we assume it is true, then we know that the delta 
function in eq. (6.7) tells us that 

r(.") oc^^(r-f) , (6.31) 
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where f stands for the maximum of the time parameters of all internal line propagators. The 
proportionality factor above does not depend on r. Thus, the only way that r appears on 
the right hand side above is as an upper bound on the integrals over the time parameters of 
the internal lines. Let us denote the maximum of the time parameters of the external legs 
of r^""^ as 77max- If we can prove that r/max > t, then for r > r/max, the right hand side of 
eq. (6.31) will remain constant, thus proving eq. (6.6). 
So, let us prove that 

7?max > T . (6.32) 

Let us assume to the contrary: there exists an internal propagator with a time parameter 
which succeeds all external leg time parameters. That time parameter can be on either an 
internal wiggly line or on an internal straight line. If it is on a straight line, then according 
to (6.21) and (6.23), there must exist a wiggly line in the diagram which comes later. It 
cannot be an external leg, or otherwise our assumption will be violated. So, it must be that 
the latest internal time parameter must be on an internal wiggly line. However, that wiggly 
line can attach only to the internal propagator ^t^b due to eq. (6.22). And therefore, we see 
that there exists a straight internal line which has the latest internal time parameter. But 
we already saw that this is impossible. Thus, our assumption must be wrong, and therefore 
^max > ^- Therefore, eq. (6.6) holds. 

(n) 

One last property of T)- deserves our attention. For r < r^max we can prove using an 
analysis similar to the above (see Appendix B of [34]) that 

r(") = rj';) , for r < ??^ax . (6.33) 

Thus, the flow of r^"^ turns out to be rather trivial under the sharp temporal cutoff: = 
rj,"^ = S^") for r < ?/max, and jumps to T^^ = T^") for r > 7?max- 

(n) 

Now let us show that Wr depends on the statistical properties of the system only 

(n) 

at previous times. To see this, note that each IPX vertex entering in the cumulants, Wr , 
must have an external wiggly line, which corresponds to the latest time parameter (as shown 
above) entering in that IPI diagram. We must chain at least one IPI diagrams with internal 
propagators, satisfying (6.22). Thus, we end up with an external wiggly line which has the 
largest time parameter, r]x, among the time parameters in the chain. The only external 
propagator that can attach to that line is h^-ax, which has i]a > r].x, with 7]a being the time 

(n) (n) 

of the external leg of W r . Thus, Wr depends only on the past properties of the system, 
which can be seen in e.g. eq. (6.17). 

This concludes our proof that under the sharp temporal cutoff, the structure of the 
NPRG flow and its solutions, the Helmholtz hierarchy, is consistent with causality; and that 
the auxiliary field, X) is unobservable. This property is preserved under truncation of the 
hierarchy as long as any truncation to the NPRG (or Helmholtz) hierarchy is consistent with 
all causality conditions discussed above. 

6.5 Summary and discussion 

In this section we integrated the NPRG flow equations for the CDM using the sharp tem- 
poral cutoff, proposed in [21] for dealing with non-equilibrium quantum field dynamics. The 

■^^Eq. (6.6) was proved in [34] using similar arguments. 



- 32 - 



solution is what we call the Helmholtz hierarchy, the first three orders of which are given by 
equations (6.14,6.15,6.16) supplemented by (6.9,6.10,6.18,6.27). 

(n) 

We showed that the presence of the cutoff r in , implies that the time parameters in 
all internal loops are bounded from above by r. If r = ry/, we directly obtain the anticipated 
result: ri%, = 5("). Thus, by moving r we interpolate between the bare action S for the 
CDM and the fully integrated effective action, T, by including gravitational interactions up 
until the cutoff in time, r. The flow of r^"^ turns out to be a simple jump from 5*-"-* to r*-"-* 
at T = 7?max (the largest time parameter of the external legs of F^")). 

Truncating the HH at ©(F^")) for n > 4 is in general a better approximation than 
truncating it at 0(e") since F^") is a sum of terms which are O ((if/^)") ~ 0{e"') or higher. 
However, a simple F^") truncation is not physical, as there is no manifest physical ordering 
parameter. So, the way one should approach the HH is by first truncating at ©(F^")) for 
n > 4, and then truncate at ©(e"). For n < 4, one should be careful with the bare values 
contributing to F^") (e.g. at lowest order F.^jj = T^^.^jj ~ Ois)). 

We showed explicitly that the truncated NPRG and Helmholtz hierarchies obey causal- 
ity, equations (4.7) and (6.4), to all orders. Moreover, we showed that the generated auxiliary 
field is unobservable, as required by eq. (4.8). 

Combining eq.s (6.9,6.14,6.16) we find the equation of motion for /: 

Lax fx LCaxyfxfy l^a — K^xyGxy i (6.34) 

where we used that the right hand side of eq. (6.16) vanishes exactly as per eq. (A.l) in 
Appendix A. Thus, the above equation is exact, and note that it is identical to the first 
extended BBGKY equation, eq. (3.2). 

The two-point correlation function can also be solved for easily to obtain (using (4.13) 
and (6.17)): 



T 




T 



^0\T^y-f) + (terms including r/^^,rxxx) (6.35) 

where the first equality defines the black blob as the amputated two-point correlation func- 
tion. The diagrams above are drawn so that time flows along the horizontal axis in each 
diagram. That is done to highlight the causal structure of the diagrams. Thus we see that 
the two-point function is generated by two distinct contributions: The first term in (6.35) 
corresponds to the initial conditions propagated forward in time; while the rest of the terms 
show the generation of power through mode coupling. A result with the same structure was 
obtained in Renormalized Perturbation Theory (RPT) (e.g. eq. (9) in [35]). As an example 
of how mode coupling generates power in the two-point function, in the second term above we 
can see that two pairs of modes are generated (at the black blobs) either by propagating the 
initial conditions or by nonlinearities, and then those pairs of modes interact gravitationally 
(at the gray vertices) to dump power to the two point function. 
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The above equation can be simplified using the 9 function in eq. (6.35) as follows: 

Gab = ■^iiax^xyRyb (6.36) 

= 47^^12^13^24 ( f f y ) ^(ri2 - T34) + {x^y) + 



'-xy = '*-'v^l2>-^13>-^24 2^ riaJs^Xy 

+(terms including ^ j^^ and ^-xxx) ■ (6.37) 

As in (6.35) the term ^ax^xy^yb tells us how the initial 2-point function is propagated 
forward in time by the response function, R. The term Rax^xyRyb tells us how much 
power is generated from mode coupling due to the nonlinear evolution of CDM. Thus, IT is 
usually referred to as the "self-energy" of the system. Note that not surprisingly eq. (6.36) 
has an identical structure to the equation for the two-point function in RPT (the second 
equation in Fig. 10 of [8]). 

Truncating the above equation at 0{e^), the first brackets in (6.37) become i^j/34, and 
the last line of that equation drops out. Thus, one obtains 

^xy = 4:Kxi2GisG24,Ky3/i . (6.38) 

Equations (6.36) and (6.38) are equivalent to eq. 45 of [19] (after using his eq. 47; keeping the 
contribution of the initial conditions; and changing to his set of variables). One important 
difference from [19], however, is that to be consistent, we have to truncate (6.36) and (6.38) 
at 0{e^). 

The linear response function (referred to as the "nonlinear propagator" in [8]) is given 

by 

xd^{rax-f) + (terms including rj^^,r^^^) ^^g gg^ 

Note that the label a in the diagram above lies on a gray wiggly line. The K term on the 
left hand side shows us how the disturbance generating the linear response couples to the 
average CDM distribution, /. The term represented by the diagram, shows how the response 
function is modified by mode coupling, where a pair of modes (generated at the black blob) 
interact with the response function at the vertices (gray vertices). The structure of the above 
equation is the same as that in RPT (the first equation in Fig. 10 of [8]). Using the causal 
structure of the vertices in analogy as we did in Appendix A, the above equation combined 
with the equation for the other non-diagonal element of VF^^^ (see (4.13)) can be written as 

{Lax — "^Kaxyfy) ^xb = ^ab + ^ax'^xb (6.40) 
{Lxa — "^Kxyafy) Rxb = ^ab + ^axRxb (6-41) 




^) ^-^34 ;/<^xi/2 -^^24-^^634 + (terms includingF.j^^ and F; 



(6.42) 

The term containing ^ ~ the "damping self-energy" - serves a similar role as IT, i.e. it 
redistributes power due to the nonlinear interactions. Note that the above expression for the 
damping self-energy is equivalent to the expression obtained by [19] (cf. his eq. 56). 
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At 0{e^), the above set of equations is identical to the result in [19] (his eq. (46), (47)), 
reducing the damping self-energy to: 



Sab = 4:Ki2aRl3G24:Kb3A ■ 



(6.43) 



Again note that to be consistent with our physical expansion, we have to truncate (6.40), 
(6.41) and (6.43) at 0{£^), unlike [19]. 

As a check of our results, let us see how one can recover the final result by Valageas 
[19] from the HH. If one does not truncate our equations with respect to e, but instead 
work to 0{P1), one recovers the result (eq. (42)-(50)) obtained by [19] using the steepest- 
descent method applied to a large-A^ expansion of the action. Thus, if instead of using 
e, we use Pl as our ordering parameter, we can recover the result of second order standard 
perturbation theory [7], which result was re-derived by Valageas [19] using his equations (42)- 
(50) and initial conditions from which important phase-space information has been dropped 
(see discussion in Section 2.2) by the Pl expansion. 

Now let us see how the HH recovers information about stream crossing. The easiest 
way to show that is to work to zero order in e, which is equivalent to the ZA. In that case 
we have 



Note that the last equation is equivalent to (2.21) as must be the case. The only difference 
between the exact two-point function (2.14) in the ZA and {fafb)c entering in A^;, is the time 
arguments. That difference is removed by the Dirac delta functions in h,ax in the expression 
for Gab above. Thus, Gab above readily reduces to (2.14). Therefore, we can see that the HH 
at zero order in e directly recovers the exact 2-point function in ZA. Thus, we can conclude 
that the HH preserves the information about stream crossing. 

From (6.44,6.45) one can also see that the high k decay kernel in the two-point function 
(2.15) is entirely due to the presence of the decay kernel in the initial conditions, Aab, where 
the decay in w is translated to a decay in k by the delta functions in 1?^. The structure of 
h,ax is such that the decay at high w of the zero order /, (2.18), and Aab will be converted to 
a decay at high k of the CDM power spectrum to all orders in e. Thus, the behavior of the 
density response function at high k anticipated in RPT [3] can be automatically reproduced 
if desired by using the exact initial conditions in phase-space. 

Next, let us comment on the relation between the extended BBGKY and the Helmholtz 
hierarchies. Combining the extended BBGKY equation for the linear response, (3.6), with 
(4.9), (4.12), and (6.29) one can recover the second Helmholtz equation, eq. (6.40) above. 
From the extended BBGKY equation for the 2-pt function, (3.3), using the same procedure, 
after a bit of algebra one can recover (6.36). 

This relationship between the first couple of equations of the (extended) BBGKY hier- 
archy and the HH should not be surprising. The (extended) BBGKY hierarchy is simply the 
hierarchy governing the evolution of W^"'\ while the HH governs the evolution of F^"). As W 
and F are related by a Legendre transformation, the two hierarchies must lead to identical 
dynamics. 

Therefore, it should be possible to extract a closure relation for the extended BBGKY 
hierarchy starting from the truncated HH. Using (6.29), one can start from the 3-vertex, F^'^^, 




Wa-Wb+ {D{ria) - D{r]b))k, 



■a 



(6.44) 



(6.45) 
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to obtain the three-point correlation function {fafbfc)c of the CDM, and the two-point hnear 
response, S{fafb)c/SCc- Dropping the 0{e^) terms which include diagrams with vertices ^^^f 
and ^xxx^ ^"P^ function we obtain: 




{fafbfc)c = +3 sym. terms 



(6.46) 



where we made explicit the value of r for the vertices F^^^ . The 2-pt linear response functions 
is given by: 




^IMkk = V + (a^b) 

^ (6.47) 

If one is to re-include the diagrams with vertices T^^j and ^xxx above two equations, 

the equations become exact. 

The interpretation of the closure relations above is straightforward. The closure rela- 
tion for the 3-pt function, (6.46), tells us that a three point function is generated by the 
gravitational coupling (gray vertex) between two modes (which originate from two distinct 
pairs of modes, generated at the black blobs), the result of which is then fed into the response 
function, which propagates the effect to a third mode. Meanwhile, the closure relation for 
the 2-point linear response function, (6.47), tells us that under a non-random forcing, the 
2-pt function varies because of the coupling of the linear response to a second mode (coming 
from a pair of modes generated at the black blob). 

Therefore, the Helmholtz hierarchy collapses to the extended BBGKY hierarchy, (equa- 
tions (3.2), (3.3) and (3.6)), with closure relations given by the three-point function: 

{fafbfc )c = 2^^KxyzGybGzc + 3 sym. terms + 0{e^) 

= 2^^K^y,{fyh),{fJ,), + 3 sym. terms + ©(e^) , (6.48) 



and the two-point linear response function: 

^xyzGyb 



^^f'^f^^^ TK^K^yzGybXc + (a o 6) + 0(e3) 



= 2^K,yz{fyft,)J-^ + (ao6) + 0(e3) . (6.49) 



Note, however, that the truncation of T at 0{e'^) still generates all possible correlation and 
response functions which enter the extended BBGKY hierarchy. Thus, the BBGKY hierarchy 
is not truncated (there is no n above which (/")c is zero) under the closure relations above. 

7 Restoring the non-Gaussianities in // 

So far we assumed that // is a Gaussian random field. However, as we discussed in Section 2.2 
and as we will see below keeping the non-Gaussian part of // will turn out to be extremely 
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important if one wants the expansion of the Helmholtz hierarchy for the CDM in e to be 
self-consistent. 

Restoring the non-Gaussian part of // is rather straightforward. The starting point 
again is the ensemble average of a functional, F[f], which can be written as: 

(F[/]) = j VfiFlfUfi) = j VfiF[f] j Vye^yf' j V f[e-'yf'ip{f[) 

= I VfiF[f] j Vye'yf'{e~^y-fi) , (7.1) 

where p{fi) is the probability distribution from which // is drawn. The term in the angular 
brackets above is simply exp(Wj[—iy]), where Wj is the Gibbs free energy for // which can 
be written as: 



Wi[j] = Y.-(f?)c- (7-2) 

n=l 

Performing the same manipulations on eq. (7.1) as we did on eq. (4.1) in Section 4.1, 
we obtain 

{F[f]) = j P0F[/]e-^[^l , where (7.3) 

S[4>] = XaLahfb - XaKabcfbfc " XafJ'a " ^Xa^abXb " (7.4) 
--^XaXbXc{fl,afl,bfl,c)c - ■^XaXbXcXd{fl,afl,bfl,cfl,d)c " ■ ■ • , 

where the first line of the expression for the bare action S corresponds identically to eq. (4.5). 
The bare action above has all the necessary properties that we required in Sections 6.1 and 
6.4 in order for the HH to lead to self-consistent and causal dynamics. 

Note that now we have modified initial conditions for Tr, sich that we have non-zero 
(5"r//(5x"' for all n. This modifies the power counting in e. However, this effect is still 
perturbative and tractable. The reason for that is because in a diagram for a given n-point 
function the legs of 6^Tj/dx^ must be attached to the propagator Rab, which in turn can 
correspond to either external legs or must attach to the bare vertex K which is 0(e). If all 
Rab correspond to external legs then to zero order that term is simply the n-point function 
in the ZA. Therefore, a given n-point function is given by the corresponding n-point function 
in the ZA plus higher order corrections. This is reasonable, as we explicitly chose the ZA as 
our zero-order solution. 

Note that none of Equations (6.35)-(6.42) except eq. (6.38) are modified, as they explic- 
itly omit the ^xxx t^rm which arises from the non-Gaussian contributions to //. Truncating 
those equations at third order we find that only the self-energy receives a modification from 
non-Gaussianities in // as follows: 



Hxy — 4:Kxl2GisG24KyS4 + 4:{fl^afl,xfl,y)cRxuRyvKzuvRzb ■ (7.5) 

Working to zero order in e we again recover the ZA. At first order we obtain 

iS(^) = iC^'^ {2K,.yff^ - lW) , (7.6) 

G^J = 2^(0) (2Kzxyff^ - L«) + 2{faMy)^^)Kz^ylCb^^) , (7.7) 
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where the expression for G^^j must be symmetrized. Here a subscript in parenthesis denotes 
the order in the e-expansion to which we work in. The first order quantity, L^^\ is the term in 
L which contains w • dw The last term in the expression for G^^^ is due to the non-Gaussian 
contributions from //. Note that ah zero-order quantities are evaluated in the ZA. 

We can also write the closure relations for the BBGKY hierarchy including the non- 
Gaussian contributions from //. Both eq. (6.47) and (6.49) remain unchanged, except they 
now have corrections of order 0{e'^). We recover eq. (6.48) but with the addition of the initial 
3-point function which can be evaluated in the ZA: 

{fafbfc )c = 'iiax^iby^icz{fl,xfl,yfl,z)c + (7.8) 

+ (2^^KxyzGybGzc + 3 sym. terms^ + 0{e'^ 

_ Sfa Sfb Sfc 
^Cx Ky Kz 



-{fl,xfl,yfl,z)c + 



+ ( 2^K,y,{fyfb)c{fzfc)c + 3 sym. terms^ + ©(e^) 



Similarly, one needs to modify eq. (6.46). Note that when evaluated to zero order in e the 
above equation reduces to the Zel'dovich 3-pt function. 

8 Particle interpretation 

The solutions for the 2-pt function provided in the previous section for different orders in e 
are numerically hard to solve, because of the phase-space n-pt functions in the ZA. 

Therefore, in this section we will rewrite our results in terms of effective particle tra- 
jectories. Those would allow one to run N-body simulations which will recover the different 
orders in e. In order to proceed, note that one can write the exact VP equation as: 



fx — ^az^ K^xyfxfy ^^x fx 



i.l) 



Using a Born-like approximation one can solve the above equation iteratively by starting with 
a solution which is close to the true solution for /. Choosing the initial / for the Born-like 
approximation to be / in the ZA, one can obtain a solution for /. That solution should 
have identical statistics as the ones obtained through the Helmholtz hierarchy, which uses 
the same expansion parameter as the Born-like approximation in the equation above. 

Indeed at first order eq. (8.1) reduces identically to the first order result in the HH, 
eq. (7.7). After some algebra at first order the above equation yields 

/'■)(., .,,) = /<i»,/d,/i,^l!!g^x (8.2) 



X (^d^,<^'^^\x',i^')+D's{q)^ ■d^ + {D-D')(^d^,^^^\x',r^')+D's{q)^ ■ 
6D{x-q- Ds{q)) 5d {v - s{q)) , 
where x' = q + D's{q), and primed quantities are evaluated at r]'. And $ is defined as 

V2^> EE 5 . 
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Note that to first order in 6, we have = —Ds, and therefore f^^^ above contains contri- 
butions which are only second order and higher in the density perturbations. Moreover, this 
ensures that the time integral above converges. 

If we chose to not restore the non-Gaussian contributions to // in the first-order equation 
in the HH, eq. (7.7), then one can check that this is equivalent to using an inconsistent ^ 
which contains a zero order piece in the overdensity. This in turn would make the time- 
integral above divergent, and the HH would have broken down even at first order. Thus, 
keeping the non-Gaussian contributions to // is crucial for the self-consistency of the HH. 

We can extract an equation of motion for particles which are described by the first order 
one-particle distribution above. To do that let us write their trajectories as 

x{q,ri) = q + Ds{q) +y{q,vi) . (8.3) 

The corresponding one-particle distribution function has the form: 

f{x,v,rj) = J d'q dD{x-q- Ds{q) - y{q, rj)) 5d - s{q) - . (8.4) 

Expanding the above expression to first order in y (denoted by a superscript (ly)) we obtain 

f^'y\x,v,r^) = -J d\ (^yiq,^).dx + ^^-d^^ (8.5) 

5D{x-q- Ds{q)) 6d {v - s{q)) . 

Comparing f^^y^ in the above equation with /^^^ from eq. (8.2) we can read off both y{q,r]) 
and if{q,r]), which are consistent with one another. So, indeed the HH hierarchy at first 
order can be treated by performing an N-body simulation using particles with trajectories^^: 

V 

xiq,^) = q + Dsiq) - J dr,' {D - D')^, '^^''^^^j^'^ {dx'^^'\x' ,r,') + D' s{q)^ , 



(8.6) 

where x' = q + D' s{q). This x{q,r]) is a solution to the following equation of motion 

D D dln(aD) f ^ ,tn\, ^ ^ , ^\ , -, / 

X -X = ^ dx^^^'(x, T]) + Dsiq) , or equivalently: (8.7) 

D D drj \ J 

x + nx + dx^^''> = -a-'^ ]'dr,'y-^^^{dx>'^^'\x',^') + D's{q)^ , (8.8) 



where we re-expressed $ with the Newtonian gravitational potential (p on the left hand side 
of eq. (8.8). One can compare eq. (8.8) with the true equation of motion given by 

x + 'Hx + dx(t) = Q . (8.9) 



^^Note that when we extracted the relevant particle trajectories, we assumed an expansion in the displace- 
ment, y, with respect to the particle trajectories in the ZA. Thus, eq. (8.6) will in fact lead to a one particle 
distribution function which features additional resummations in y with respect to the first order in HH. 
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Thus, the right hand side of eq.(8.8), which is first order in e and second order and higher in 
the overdensity, corrects for the fact that the potential <j) on the left hand side is evaluated 
in the Zel'dovich approximation. 

The above equation of motion, eq. (8.7) or (8.8), is not very convenient for implementing 
in an N-body simulation, both because of the non-standard velocity it uses (cf. eq. (8.7)), 
and because generalizing it to higher orders is not straightforward. One can modify a bit 
the expansion scheme above to remove both of these problems. This is simply achieved by 
replacing {D — D') in eq. (8.6) by J"^ dr]"!/ aij]") and dropping the logarithm: 

x{q,i^) = q + Ds{q)- (8.10) 


One can check that when evaluated with the full and with the exact x' (i.e. not with 
x' evaluated in the ZA), the above equation is the exact solution to the true equation of 
motion, eq. (8.9). This solution is written in a form which contains explicitly the ZA, and is 
manifestly rather close to the e-expansion^^ . 

Writing the equation of motion eq. (8.9) and its solution eq. (8.10) to higher orders in 
(V<^ + Ds) is straightforward: 

+ nx'^'' = -5a,(^^("-^)(a;^"^) , (8.11) 

where a subscript (^ n) implies that the quantity is a sum of all orders up to and including 
n. Above one has a choice of m = n or m = n — 1. If we are to stick as close as possible to 
the HH, then we should choose m = n — 1. However, m = n may offer faster convergence - 
this is something which needs to be investigated further numerically. The zero order is given 
by (^(°) and £c(°) in the ZA. 

Equation (8.11) represents an iterative scheme for improving the Zel'dovich approxima- 
tion. It is inspired by the e expansion, and if required can be modified to conform exactly to 
the £ expansion in the same way in which we obtained the first order solution, eq. (8.6). 

The above iterative scheme for solving the VP equation is well suited for implementing 
in N-body codes. This can be achieved in the way outlined in Section 1. 

9 Summary 

In this paper we derived a self-consistent hierarchy (the Helmholtz hierarchy) of partial 
differential equations, governing the evolution of the phase-space statistics of cold dark matter 
in the absence of baryons. The Helmholtz hierarchy (HH) has a physical ordering parameter, 
which interpolates between Zel'dovich dynamics and fully-hedged gravitational dynamics. It 
is schematically given by the fractional difference between the acceleration of test particles 
as given by a Zel'dovich-type approximation (i.e. the acceleration is assumed parallel to the 
velocity at an intermediate moment in time), and their corresponding true acceleration due 
to gravity. 

We showed that the HH preserves information about stream crossing and we argued 
that it automatically generates the decay expected at high k of the density response function 

■^^Indeed, starting from the VP equation and expanding in our variable y one can show that the two are 
related by resumming certain higher order terms in e. 
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[3]. However, unlike RPT and analogous approaches, the resulting non-linear power spectrum 
is expanded around the ZA, possibly allowing us to construct better models for the mildly 
non-linear (and possibly non-linear) regime. 

Under a sharp truncation of the HH all n-point correlation functions of CDM are gen- 
erated, in contrast to the BBGKY hierarchy. Combining this result with the presence of 
a physical ordering parameter, we find that the HH ameliorates the closure problem of the 
BBGKY hierarchy. 

We derived the HH for CDM using the functional renormalization group with a temporal 
cutoff [21], so that causality is built-in from the start. We proved that causality is indeed 
preserved to all orders — a nontrivial result^^, which is important since structure formation 
is an inherently out-of-equilibrium process. 

The HH has several advantages over performing a numerical N-body simulation and 
then taking spatial averages to reproduce the ensemble averages. First, the correlation and 
response functions are smooth (over the domain where they are nonzero), unlike the one- 
particle distribution function, or in the case of N-body simulations, the density field. Second, 
the effect of the short-scale modes on the long-scale modes (and vice versa) is readily taken 
into account by the HH. Third, the HH may offer the opportunity to develop better analyt- 
ical or semi-analytical templates for fitting the nonlinear part of the CDM power spectrum 
obtained from numerical simulations. 

The HH hierarchy is easy to solve numerically at zero order, since it corresponds identi- 
cally to the ZA. However, going to higher orders requires a good numerical handle of the full 
phase-space n-point functions in the Zel'dovich approximation. To go around this problem, 
we proposed an iterative scheme which closely follows the e-expansion. The scheme uses 
successive N-body simulations which improve upon the Zel'dovich approximation. However, 
we postpone such a numerical investigation to a future study. 

A Simplifying the contribution of r*^"+^) to F*^") 

In the solutions to the NPRG hierarchy with a temporal cutoff, the time parameters in all 
IPI diagrams must come earlier than the external legs of the corresponding IPI diagram. 
This comes as a direct consequence of the delta function in (6.7). Examples of this can be 
seen in eq.s (6.15,6.16). 

Thus, the following equation holds for effective actions satisfying (6.21), (6.23) and 
(6.25): 




(A.l) 



where r is the maximum of the time parameters running inside F^"); dashed lines represent 
either wiggly or straight lines. To derive the above equation we used the causality conditions 

■^^The usual calculations in statistical mechanics assume that the system under study is in equilibrium. 
In such cases, the resulting equations of motion usually have either no dependence on a time parameter, or 
causality is not manifest. 
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(6.23) combined with (6.22), and the 0-function restrictions written explicitly in the Feynman 
diagram above. Up to a proportionality constant, the left hand side is the general form of 
the r("=™-+2) term in the equation for T^"^\ Therefore, using tI^-'^^ = for the CDM, from 
(A.l) we find that the Pr""*"^^ term vanishes in the expression for r^"^ for 7i > 1. 

Acknowledgments 

The author would like to thank Matias Zaldarriaga for many helpful conversations. 



References 

[1] S. Cole, S. Hatton, D. H. Weinberg and C. S. Frenk, "Mock 2dF and SDSS galaxy redshift 
surveys," Mon. Not. Roy Astron. See. 300, 945C (1998) [arXiv:astro-pli/9801250]. 

[2] R. Angulo, C. M. Baugh, C. S. Frenk and C. G. Lacey, "The detectability of baryonic acoustic 
oscillations in future galaxy surveys," Mon. Not. Roy. Astron. Soc. 383, 755 (2008) 
[arXiv:astro-pli/0702543]. 

[3] M. Crocce and R. Scoccimarro, "Memory of Initial Conditions in Gravitational Clustering," 
Phys. Rev. D 73, 063520 (2006) [arXiv:astro-pli/0509419]. 

[4] Ya. B. Zel'Dovich, "Gravitational instability: An approximate theory for large density 
perturbations." Astronomy & Astrophysics 5, 84 (1970) 

[5] P. Valageas, "Using the Zeldovich dynamics to test expansion schemes," arXiv:0706.2593 
[astro-ph] . 

[6] D. Baumann, A. Nicolis, L. Senatore and M. Zaldarriaga, "Cosmological Non-Linearities as an 
Effective Fluid," arXiv:1004.2488 [astro-ph.CO]. 

[7] B. Jain and E. Bertschinger, "Second Order Power Spectrum and Nonlinear Evolution at High 
Redshift," Astrophys. J. 431, 495 (1994) [arXiv:astro-ph/9311070] 

[8] M. Crocce and R. Scoccimarro, "Renormalized Cosmological Perturbation Theory," Phys. Rev. 
D 73, 063519 (2006) [arXiv:astro-ph/0509418]. 

[9] P. Valageas, "Large-N expansions applied to gravitational clustering," arXiv:astro-ph/0611849. 

[10] S. Matarresc and M. Pietroni, "Baryonic acoustic oscillations via the renormalization group." 
Mod. Phys. Lett. A 23, 25 (2008) [arXiv:astro-ph/0702653]. 

[11] N. Afshordi, "How well can (renormalized) perturbation theory predict dark matter clustering 
properties?," Phys. Rev. D 75, 021302 (2007) [arXiv:astro-ph/0610336]. 

[12] P. Valageas, "Impact of shell crossing and scope of perturbative approaches in real and redshift 
space," Astron. Astrophys. 526, A67 (2011) [arXiv:1009.0106 [astro-ph.CO]]. 

[13] J. Carlson, M. White and N. Padmanabhan, "A critical look at cosmological perturbation 
theory techniques," Phys. Rev. D 80, 043531 (2009) [arXiv:0905.0479 [astro-ph.CO]]. 

[14] P. McDonald, "How to generate a significant effective temperature for cold dark matter, from 
first principles," arXiv:0910.1002 [astro-ph.CO]. 

[15] S. Pueblas and R. Scoccimarro, "Generation of Vorticity and Velocity Dispersion by Orbit 
Crossing," Phys. Rev. D 80, 043504 (2009) [arXiv:0809.4606 [astro-ph]]. 

[16] T. Matsubara, "Rcsumming Cosmological Perturbations via the Lagrangian Picture: One-loop 
Resuhs in Real Space and in Redshift Space," Phys. Rev. D 77, 063530 (2008) 
[arXiv:0711.2521 [astro-ph]]. 

[17] P. Catelan, "Lagrangian Dynamics In Nonflat Universes And Nonlinear Gravitational 
Evolution," Mon. Not. Roy Astron. Soc. 276, 115 (1995) [arXiv:astro-ph/9406016]. 



- 42 - 



p. J. E. Peebles, "The large-scale structure of the universe," Princeton University Press, (1980). 

P. Valageas, "A new approach to gravitational clustering: a path-integral formalism and 
large-N expansions," Astron. Astrophys. 421, 23 (2004) [arXiv:astro-ph/0307008]. 

M. Croccc, S. Pueblas and R. Scoccimarro, "Transients from Initial Conditions in Cosmological 
Simulations," Mon. Not. Roy. Astron. Soc. 373, 369 (2006) [arXiv:astro-ph/0606505]. 

T. Gasenzer and J. M. Pawlowski, "Functional renormalisation group approach to 
far-from-equilibrium quantum field dynamics," Phys. Lett. B 670, 135 (2008) [arXiv:0710.4627] 

S. Bharadwaj, "The Evolution of Correlation Functions in the Zel'dovich Approximation and 
its Implications for the Validity of Perturbation Theory," Astrophys. J. 472, 1 (1996) 
[arXiv:astro-ph/9606 12 1] 

C. P. Ma and E. Bertschinger, "A Cosmological Kinetic Theory for the Evolution of Cold Dark 
Matter Halos with Substructure: Quasi-Linear Theory," Astrophys. J. 612, 28 (2004) 
[arXiv:astro-ph/0311049]. 

M. Davis and P. J. E. Peebles, "On the integration of the BBGKY equations for the 
development of strongly nonlinear clustering in an expanding universe," Astrophys. J. Suppl. 
Series 34, 425 (1977). 

P. C. Martin, E. D. Siggia and H. A. Rose, "Statistical Dynamics of Classical Systems," , Phys. 
Rev. A, 8, 423 (1973). 

R. Phythian, "The functional formalism of classical statistical dynamics," Journal of Physics 
A: Mathematical General, 10, 777 (1977). 

A. Kamcncv, "Many-body theory of non-equilibrium systems," arXiv:cond-mat/0412296. 

M. Kardar, "Statistical Physics of Particles," Cambridge University Press, (2007). 

M. E. Peskin and D. V. Schroeder, "An Introduction To Quantum Field Theory," 
Addison-Wesley, (1995). 

C. Wetterich, "Exact evolution equation for the effective potential," Phys. Lett. B 301, 90 
(1993). 

J. Bergcs, N. Tetradis and C. Wetterich, "Non-perturbativc renormalization flow in quantum 
field theory and statistical physics," Phys. Rcpt. 363, 223 (2002) [arXiv:hcp-ph/0005122]. 

K. G. Wilson and J. B. Kogut, "The Renormalization group and the cpsilon expansion," Phys. 
Rcpt. 12, 75 (1974). 

J. M. Pawlowski, "Aspects of the functional renormalisation group," Annals Phys. 322, 2831 
(2007) [arXiv:hep-th/0512261]. 

T. Gasenzer, S. Kefiler and J. M. Pawlowski, "Far-from-equilibrium quantum many-body 
dynamics," European Physical Journal C, 223 (2010) [arXiv:1003.4163] 

M. Crocce and R. Scoccimarro, "Nonlinear Evolution of Baryon Acoustic Oscillations," Phys. 
Rev. D 77, 023533 (2008) [arXiv:0704.2783 [astro-ph]]. 



-43- 



